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The technique of density matrix equation (DME) for a small system interacting with a bath is 
explained in detail. Special attention is given to the nonsecular DME that is needed in the vicinity 
of overdamped tunnelling resonances in molecular magnets (MM). The relaxation terms of the 
DME for MM are represented in the universal form that does not employ any unknown spin-lattice 
coupling constants and absorbs the information about the spin Hamiltonian in the exact basis states 
and transition frequencies. This makes adding new types of anisotropy easy and error-free. The 
Mathematica code is available from the author. 
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I. GENERAL THEORY 

Density matrix is used to describe properties of a system that is a part of a larger system with which it interacts. 
Whereas isolated systems (such as the above mentioned larger system) can be described by a Schrodinger equation, 
systems which interact with their environments cannot. Starting from the Schrodinger equation for the isolated whole 
system, small system + environment (or bath), and eliminating the environmental variables, one can, in principle, 
construct an object that can be used to calculate observables of the small system in a short way. This object is the 
density matrix of the small system. Of course, integrating or taking matrix elements in two steps, at first over the 
environment and then over the small system is not a big simplification. This approach becomes really useful if one 
obtains a closed equation of motion for the density matrix of the small system, the density matrix equation (DME). 
This is possible if the interaction between the small system and its environment is small and can be considered as 
a perturbation, and the small system does not strongly perturb the state of the environment. The derivation of the 
DME for a bathed small system will be presented in Sec. [Til Here the necessary components of the formalism will be 
introduced. 



A. From the wave function to the density matrix 



The general wave function or state \ip) of an isolated quantum system can be expanded over a set of complete basis 
states \^ m ) as 

|*)=2cm|*m). (1) 

m 

The coefficients c m completely characterize the state |^) and can be used to calculate physical quantities A described 
by corresponding operators A: 

A=(A) = = ]T C m <(*„|i|* m ) = C ™ C n A nm- (2) 

mn mn 

One can define the density matrix p corresponding to quantum-mechanical states (pure states) of our small system 
by its matrix elements as 

{P}mn = Pmn = c mC n: (3) 

then Eq. ((2]) becomes 

A = ^P mn A nm . (4) 



■nil I 



The density matrix satisfies the normalization condition 



Tr{p} = V Pmm = l. (5) 



This condition follows from Eq. ([3]) for the pure states but it holds in general, too, since the average of the unity 
operator A mn — S rnn should be 1. Additionally, the condition 

Ei^«i 2 = 1 ( 6 ) 

mn 

is satisfied for the density partix of pure states, Eq. 

For systems interacting with their environment, observables A still are given by Eq. ([4]), although the coefficients 
p nm in general do not reduce to products as in Eq. ([3]). For the whole system including the environment, one can use 
basis states that are direct products of those of the small system \ip m ) and those of the environment |</> ro ): 

|*m«> = IVO ® \4>J = \^m<f>J- (7) 
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The quantum mechanical states of the whole system (considered as isolated) can be written, similarly to Eq. ([1}, m 
the form 

I*) =^2C mza \V mza ). (8) 

tnvj 

The expression for the observable A of the small system becomes 

A = (V\A\*) = E C* mv ,C nv ,,{* mw \A\* nza ,) = E C; n ^ m ,^JA\i, n ){^\^,). (9) 

For the orthonormal set of |</> ro ) this can be rewritten in the form of Eq. ((4]) with p mn p^„, where (f mn is the 
reduced density matrix of the subsystem given by 

Pmn — ^ t Cm^Crm = ^ t Pm-ca.nzo- (1^) 
zu w 

Since, in general, C mro do not split into the factors depending on m and w, the reduced DM p s mn is not a product of 
its wave-function coefficients, p s mn ^ c m c* . Obviously p s mn in Eq. (flU)) depends on the state of the whole system. One 
can check that it satisfies the normalization, Eq. |5j). On the other hand, Eq. (J6j) is not satisfied, in general. 



B. Entangled states and quantum statistics 



It would be wrong to think that when the interaction between the small system and the environment becomes very 
weak, the density matrix of the small system, Eq. (fTOjl . simplifies to Eq. ([3]). If the coupling is vanishingly weak, 
quantum states of both the small system and the environment are well defined. Thus both the small system and the 
environment could be in their pure quantum mechanical states. However, there are pure states of the whole system 
that do not consist of pure states of the two subsystems. Such states can be called entangled. For instance, if both 
subsystems are two-level systems with the states l^i), ^2) an( i l^i)) 1^2)1 respectively, and the whole system is in the 
state 



I*) =C 11 |V 1 1 )+C22|^2^2> (11) 

with |Cn| 2 + I C22 1 2 = 1, Eq. (fT0| with C mro = C mm S mm yields the diagonal density matrix 

Pmn = |Cmm| &mn- (12) 

This density matrix satisfies Eq. ([5|). On the other hand, diagonal p B mn does not correspond to any pure state, cf. Eq. 
©. Thus Eq. © is not satisfied, 

EKJ^ElO™! 4 ^ 1 - (is) 

mn m 

In particular, for |Cn| 2 = | C22 1 2 = 1/2 the result is 1/2 instead of 1. This is unrelated to the interaction between the 
two subsystems and it alone mandates using the density matrix rather than the wave function for a subsystem in the 
general case. 

In a non-entangled state, in the absence of interaction, the wave function of the whole system factorizes: 

1*) = \<t>) = EcmhuE d -i^>- ( 14 ) 

m zd 

If the whole system and both subsystems are described by density matrices, the total density matrix under the above 
conditions factorizes such as 

PmzD .tlzd' PmnPzDZD 1 

(15) 

for a small system and the bath. This factorization means that both subsystems are completely independent. If they 
are prepared at the initial moment independently from each other, one cannot expect their entanglement that requires 
a special care. As the time goes, interaction between the subsystems can cause some entanglement that, however, 
should remain small if the interaction is weak. The measure of entanglement is the mismatch in Eq. (| 15[) . where p s mn 
is defined by Eq. (TIT))) and is defined by similar tracing out the variables of the small system. Note that both 
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rhs and lhs of Eq. (|15p are properly normalized according to Eq. ([5]) and its variants. As a number measuring the 
entanglement, one can use, e.g., 



Ent 



n 2 - 1 



^2mn \Pmn\ ^2zjtu' 


t zuzu' \ 




Pm-cu,n-DJ f 


2 



(16) 



that is related to Eq. ©. Here n is the number of states of each subsystem that is assumed to be common. If the 
whole system is in a pure state, the denominator is equal to 1. In particular, for ^ given by Eq. (jlip . n — 2, one has 
p s mn given by Eq. (O and similarly = |C roro | 2 S wzn >, so that 



4 

Ent = - 
3 



\a 



22 



(17) 



In the case \C\i\ — | C22 1 =1/2 entanglement reaches its maximal value 1. 

The coupling between subsystems also results in impossibility to describe them in terms of wave functions. The 
coupling causes slow transitions between energy levels, so that the total energy of the whole system is conserved. 
For the small system this means that it spends time in its different quantum mechanical states (in the interaction is 
weak and they are well defined), so that over a large time its state is a mixture of different states rather than a pure 
quantum mechanical state. Alternatively one can think about an ensemble of bathed small systems being in different 
pure states. An example of a mixture state is the thermal-equilibrium state 
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Em 



(18) 



where Z s = ^2 m e E ™/( k B T ) is the partition function. This formula holds if the basis states \ip m ) are eigenfunctions 
of the small-system's Hamiltonian H s , that is a natural choice. We will see below that in the course of temporal 
evolution described by the DME, the initial density matrix approaches the form above. 

In fact, the environment of the small system also can be weakly coupled to some super-bath, so that the small system 
+ environment is not an isolated system and it also should be described by the density matrix instead of the wave 
function. The role of the super-bath is just is to ensure that the whole system is not in a pure quantum state but 
in a mixture state, whereas the coupling between the bath and super-bath is assumed to be very weak and is never 
considered explicitly. In most cases the state of the environment (bath) is the thermal equilibrium that is described 
by the density matrix similar to Eq. (| lSj) . Since the bath is much larger than the small system, the weak coupling to 
the latter won't drive it out of the equilibrium. 



C. Density matrix as an operator 

The formalism of quantum mechanics allows one to consider the density matrix as an operator. More precisely, one 
can introduce the density operator 

P = ^2Pmn\^m){' i l J nl ( 19 ) 
mn 

where \tp m ) i s a complete orthogonal set of states. The density matrix consists of matrix elements of the density 
operator: 

Pmn = WUpIVOj (20) 

so that both DM and DO contain the same information and are equivalent. Using the density operator allows one to 
put formulas in a more compact form without subscripts. 

In particular, the expectation value of an operator A can be obtained as a trace of Ap over any complete orthogonal 
set of states \ip n ). The calculation is especially simple if one uses the set of states in the definition of p: 

(A) = Tl{Ap} = J2(tn\Ap\lP n ) = J2^n\^ m )(^ m \m n ) =J2 A ™Pmn, (21) 

n mn mn 

and the result coincides with Eq. ([I]). It can be proven that the trace of an operator is independent of the choice of 
the basis in Tr {. . .} and operators can be cyclically permuted under the trace symbol, so that Tr <^ Ap > — Tr < pA \ . 
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If the system interacts with its environment described by the basis \4> m ), the total DO has the form 

P = Pmw.wlVwOW'T.'&u'l = Pmw.wlV'mKV'nl ® lO^w'l- ( 22 ) 

mm ,nzu' rn-co,n-cu' 

Calculating any observable A that corresponds to the small system can be done in two steps: First calculating the 
trace over the variables of the bath and then calculating the trace over the basis states of the small system using Eq. 
(f2"Tj) . The first step yields the reduced density operator for the small system 

p s = Tr b p = Tr ro p (23) 

that has the form of Eq. (fl"9)) with 

Pmn / .Pmtj.nzo' 

(24) 



If the whole system is in a pure state, this formula coincides with Eq. JTUJ). 

The density operator corresponding to the thermal equilibrium of the small subsystem in contact with the bath has 
the beautiful form 

(-&)■ (25 > 

Because of applications in quantum statistics density operator is also called statistical operator. Defining the density 
matrix p s mn with respect to the basis of eigenstates of H in Eq. (|20|) . one arrives at Eq. (fl"8|) . The eigenstate basis is 
the most convenient for calculating thermal averages of physical quantities. However, as pointed out above, this could 
be done with the help of any other basis. If the whole system consisting of the small subsystem and the bath is in 
contact with a super-bath, the equilibrium statistical operator of the whole system has the form similar to Eq. (|25|) . 

For a system in a pure state, the density operator defined by Eq. ([!§)) can be with the help of Eqs. (J3]) and §T§ 
rewritten in the form 

0=|V><# ( 26 ) 

D. Temporal evolution and interaction representation 

Temporal evolution of the DM or DO of an isolated system such as the small system + bath obeys the equation that 
follows from the Schrodinger equation. If, for instance, the whole system is in a pure state |^), its density operator is 
given by 

P=I*X*I (27) 
c.f. Eq. (|26[) . Then with the help of the Schrodinger equation and its conjugate 

»*^|*> =#|*>, -ih^\ = (*\H (28) 



one obtains the quantum Liouville equation 



H,p ■ (29) 



It should be noted that this equation is not an equation of motion for an operator in the Heisenberg representation 
that has another sign. The DO consists of states that have their own time dependence, unlike Heisenberg operators 
whose time dependence is borrowed from the states. In the presence of the super-bath the system is not in a pure 
state but rather in a mixture of different states I*), 



P 



= 5>*I*X*I- (30) 



Fortunately, Eq. (|29[) is valid also in this case. Remember that Eq. (|29[) describes the whole system, small system + 
bath, although we dropped the index "tot" for brevity. Our final aim is, however, to obtain a density matrix equation 
(DME) for the density matrix of the small system alone with the bath degrees of freedom eliminated. 
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Let us break up the Hamiltonian of the system into a "simple" Hamiltonian Hq and perturbation or interaction V: 

H = H a + V. (31) 

In the absence of V, evolution of quantum states \^>) is described by the unitary evolution operator Uo(t): 

|*>t = Uo(W) , (¥| t = WoUfo), (32) 

where |^)o corresponds to the starting moment t = 0. We call Uo(t) bare evolution operator. From the Schrodinger 
equation (|28|) the equations of motion for Uq and its conjugate follow : 



If i?o is time independent, the solution of these equations is 



U (t) 



-iHot/H 



JHot/H 



(33) 



(34) 



Of course, one can write similar formulas with the total Hamiltonian H as well, U (t) being the full evolution operator. 
The idea of using Uo(t) is to split the nontrivial part of the evolution due to the perturbation V from the trivial 
evolution due to Hq. To effectuate this, one can introduce the density matrix in the interaction representation 



p(t)i = ul(t)p(t)u (t). 

The equation of motion for p(t)i follows from Eqs. (|29|) and (|33|) : 



(35) 



~dt 



ih 



H,p\ U a {t) + ul{t)p{t)H U 



dt 



= -U^H p(t)U (t) + U^(t) 



= W) 



v,p 



U (t). 



Inserting f/o(t)^o(*) = 1 [ see Eq. dHS> ] between the operators and defining 



this equation can be brought into the form 



Vi{t) = ul(t)VU (t), 



i ri- 



ot 



Vi,Pi 



(36) 



(37) 



(38) 



One can see that the temporal evolution of the density matrix in the interaction representation is governed by the 
interaction only. This facilitates constructing the perturbation theory in V. 
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II. THE DENSITY MATRIX EQUATION 
A. From the full to reduced DOE 

In this section the equation of motion for the reduces density operator of a small system s weakly interacting with a 
bath will be obtained. We will be following the method of Karl Blum [l[ that is most practical, although not rigorous. 
Rigorous methods using the projection operator technique 2, 3] lead to the same result with much greater efforts. 

The Hamiltonian can be written in the form 



H = H + V, H = H S + H h , V = H s - h . 
Eq. (|38p for the DO of the whole system can be integrated over the time resulting in the integral equation 



Pit) i = p(o)j 



dt' 



v(t'h,p(t r h 



Inserting it back into Eq. (|38|) one obtains the integro-diffcrential equation 

d 



V{t)i,p{$)i 



dt' 



V(t)j, V{t') Il p{t') 1 



(39) 



(40) 



(41) 



which is still an exact relation. Although this equation seems to be more complicated than the initial equation (|38|) . 
it is convenient for perturbative treatment of the interaction V . If the small system is not entangled with the bath 
in the initial state, that is a natural assumption, a small interaction cannot cause a significant entanglement as well. 
Thus the density operator of the whole system nearly factorizes, see Eq. ([151) . Here we additionally use that the bath 
is at thermal equilibrium that cannot be noticeably distorted by a weak interaction with the small subsystem: 



p(t)i = h(t)ip c b q = ft>(*)-^- ex P ( ~t~^ 



(42) 



Note that this approximation cannot be done in Eq. (|40|) since it leads to disappearance of the effect of interaction. 
To properly describe this effect, one has then to account for the corrections to Eq. (|4"2"|) in the first order in V. This 
way is inconvenient. To the contrary, Eq. (|41[) already contains quadratic V terms that capture the main effect of 
interaction. Here taking into account small corrections to Eq. (|42[) can bring only irrelevant small terms. Now using 
Eq. (|42p one can transform Eq. (|41[) into the density operator equation for the small system by making a trace over 
the bath variables according to Eq. (f2"3"|) : 



d 
dt 



-TH 



v(t)i,pMiK q 



i 

¥ 



dt'Tr h 



V(t)j, V(t') l7 p s (t')jp e h q 



(43) 



For the couplings V we consider here, the first linear-F term disappears. 

The equation above is an integro-differential equation with integration over preceding times, t' < t, in the rhs. Such 
equations are called equations with memory and they are difficult to solve directly. In our case, however, the problem 
simplifies. The time evolution of p s (t)i In Eq. (|43|) is slow since it is governed by the weak interaction between the 
system and the bath. On the other hand, t' dependences of the other terms in the integrand (the kernel) are governed 
by H an d thus they are fast at the scale of p s {t)i- The analysis shows that the kernel in Eq. (|4"3"|) is localized in the 
region \t — t'\ < l/w max , where w ma x is the maximal frequency of the bath excitations. Thus in the integral over t 1 one 
can make the short-memory approximation p s {t')i P s {t)i after which the time integral can be calculated explicitly. 

Returning to the original reduced density operator 



p s {t)^Uo{t)p s {t)iUl{t), 
c.f. Eq. (|35[) . and computing the derivative with the help of Eq. (f3"3"| . 



H s ,p s (t) 



+ Uv(t)[^ ft p s (t)i)Uo t (t), 



one obtains 



dt 



H B ,p B (t) 



1 



dt Trb 



V. 



Vit'-t^PsWpt 



(44) 



(45) 



(46) 
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or, with t = t — t' and dropping the subscript s, p s {t) =>■ 

j t p(t) = - 1 - [H s ,p{t)\ +Rf>(t), 
where the first term with commutator is the conservative term and 



V, 



V(-T) I ,p(t)exp - 



describes the relaxation of the small system. Here, for time-independent problems, 



V{t) i = UI{t)VU q {t) 



— e iH Ty e -iH T 



(47) 



(48) 



(49) 



If the Hamiltonian of the small system H s depends on time, this time dependence is typically slow in comparizon to 
the frequency of the bath excitations w max , so that during the short times r ~ l/oj max the change of H s is negligibly 
small. Thus one can simply use Eq. (|49|) with H s = H s (t). 



B. Prom the DOE to DME 



Now one can go over from the density operator p to the density matrix p mn using Eq. (|20j) and the notations 

H s , mn = (i/J m \H s \ip n ), V mw<nvj > = (l/} m <l> m \V\lp n <f> w ,), (50) 
where \4> m ) are eigenfunctions of H^. Then the partition function of the bath becomes 



The conservative term of the resulting DME has three different forms in three different cases. If the small-system 
states \tp m } are time independent and form a so-called natural basis unrelated to H s , the DME has the form 

^[j-Pmn = {H s ,mlPln ~ Pml H s,ln) + (V'ml^lV'n)- ( 52 ) 

I 

The natural basis is inconvenient for the evaluation of the relaxation term since the relaxation of the small system 
takes place not between the states \ip m ) but between the eigenstates of H s : 



H s \x a ) 



£*\X a ), 



e-* HsT/h \x a ) = e-' l£aT/h \x c 



(53) 



The basis of \x a ) w iH be called diagonal basis since the Hamiltonian matrix (Xal-^slXa) — £ a$af3 is diagonal. In the 
diagonal basis the DME has the form 



-j t P a p = -iUccPPafi + {x a \Rp\Xp), 
where u) a p are transition frequencies between the energy levels of the small system, 

hUaf3 = £q — £/3, 



(54) 



(55) 



and the relaxation term can be conveniently evaluated (see next section). If H s depends on time, one can use the 
adiabatic basis of the states |XaM) defined as 



Hs(t)\x a (t))=e a (t)\ Xa (t)). 



In this basis p a/3 acquires additional non-adiabatic terms: 



d Pa/3 = ^(x a \f>\Xp) = (x a \f>\Xp) + (X a \p\xp) ~ ^(Xa\ 



dt 



H s ,p 



\Xp) + (xJRPlxp) 



(56) 
(57) 



(58) 



As argued at the end of the preceding section, calculation of the relaxation term (Xa\Rp\Xd) is n °t complicated by 
the time dependence of H s . This calculation will be done in the next section where Eq. (|54"|) will be used for brevity. 
The terms due to the time dependence of H s in Eq. ([55)) can be added if needed. 



C. DME in the diagonal basis 
Using Eqs. (15411 and (|48|) and inserting summation over intermediate states yields 



d 1 '* 



di p ap - -tu, afl p af} h2 j ^ 4 

a' j3' vjvj'vj" 



dr 



- E 



{{xMV\ Xa ,<t>^){xMe- iA ^/ny e ^ 
- (xMe- i6 °T/ n Ve i6 °r/ h \x a ,<t> a ,)^^ 

+ { Xc AJpe- kJ{kBT ^X a ^){x a ^e-^ (59) 



Since the states are eigenf unctions of H s and H^, this simplifies to 
± Pap = -to^-^J dr± E 

ac'P vjvj'vj" 

^ B v H 1 e V a -on,a'za' V a' zv' ,f3' ru" Pfj' f3°zct"vj 

_ i(-e ,-B„„+e +E v )r/h-E m// /(kBT) v l0lo ,5 , „V». ,,a 

^ e * (X-UJ iOl'vj' rot' /U vj 'vj" v p vj" ,pvj 

~ c ' c Paa /U roro v ol'vj' ,p vj" v p vj" ,pvj f \ uu / 

and further to 

a'/r tutu' 

/ _i(-£ a /-B„,+ E g,+£;„)T/ft„-E a ,/(fe B r) 1/ T/ 

_ e i(-^'--E„'+^+-E-)^/S e -E ro '/(fe B T)y y 

+ e^^-^^'+^WV^/^-^p^y^,^^,,^} . (61) 

Since the time kernel is sharply localized, the integration over r can be extended to the interval (0, oo) . Further, the 
relaxation of the small system that we are mainly interested in is due to the real part of the bath coupling term in 
Eq. (|61j) . Its imaginary part is in most cases only a small correction to the first (conservative) term in this equation. 
Thus taking into account the time symmetry of Fij (r) one can replace 

/ dTe<- £ °'~ E ™' +e P' +E ~) T/h => HirS (s a , - + sp + E„) (62) 
Jo 

etc. After renaming indices in Eq. (pT|) (a' A, /?' =>■ a' in the first term and j3 =>■ A, a 1 =4> in the forth term) one 
obtains the DME in the form 

= ~ iuJ a pp al3 + ^2 R a/3,a'P'P a 'l3'i ( 63 ) 
a'/3' 

where 

Raf},a'l3' = ^E j-E 6 E ™^ kBT ' ) S (s a i — £ 7 + — E m i) V aW rf m 'V'yzv' ,a'mfi fi' 
zu-gj' \ 7 

7 

+ e -^'/( fe « T ) [J ( £(3 - £/3 , + £ ro - + * (e a - e Q , + £ ro - Vw^Vp^jto} . (64) 
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Eq. (|63p can be written in a more compact form 

= ®af3., a >p>P a >p>, (65) 

where 

®a/3,a'0' = -iw a p5 aa >5j30> + R a /3,a'l3'- (66) 



In the case of time-dependent Hamiltonian one has to add the non-adiabatic terms from Eq. (|58|) to the above DME. 
The asymptotic solution of Eq. (|63[) is 

^ = i exp (n&)^ (67) 

that describes thermal equilibrium. This will be shown in Sec. Ill El 



D. Secular approximation and Fermi golden rule 

One can transform Eq. back into the interaction representation using Eqs. ([33)) and (p?U)) . In the diagonal 

basis, the relation between p a p(t) and p a p(t)i has the simple form 

P a0 (t)i = (Xje^me-^lxp) = e^p ap {t). (68) 
Computing the time derivative of p a »(t)i and using Eq. (|63[) . one arrives at the equation 

J t P a0 (t)i = £ ei^-^'^R^^p^^t)!, (69) 

where the conservative term disappeared and the relaxation term has an explicit time dependence. While the change 
of the density matrix due to the relaxation is slow, the oscillation in the terms with Lu a p ^ Lu a '0' are generally fast. 
These fast oscillating terms average out and make a negligible contribution into the dynamics of the small system. In 
general, all transition frequencies are nondegenerate, so that one can drop all terms with a ^ a' and j3 ^ f3 , if a ^ /3. 
In the equations for the diagonal terms p aa (t)i one can keep only diagonal terms with a' — . This is the secular 
approximation that tremendously simplifies the DME. In the secular approximation one has 

-^Paffi^I = <W ^ Raa, a >a>P a > a > (*) I + (1 - 5 af} ) R af 3 ,a/3 P a/3 (t) I 



8af3 ^ Ra 

a,a'a' Pa.' a' 



Simplifying Eq. (|64|) and using the Hermiticity 



OiVD ,cr VD' 



(71) 



one obtains 



and 



Raa,a' a' \ a >-£ a — ^ ^ ^ w /( )j (s a £ a > -\- E-cj E^jt ) | Vazo -uj' | — Iqq' C 7 ^) 



^ { - E e- E -/^S (e a - e 7 + E m - E m ,) \V am , JW ,\ 2 
- Yu e- E - /{kBT) 5 {ep - £ 7 + £ ro - E m ,) \V^^,\ 2 

7 

+ 2 e -E„,/(k B T) S ^ _ ^ Van^VfW,!*,} ■ ( 73 ) 
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Rearranging the terms one obtains 

Rap.aP — — r Q; 3, (74) 

where 

f a0 = T afj + - r "'« + T t 3 'P ■ ( 75 ) 

Here r Q / Q is defined by Eq. ([7^)1 and 

fa P = T7~Y. e~ E ™ /ikBT) HE^ - E„,) \V a ^, - V^,^\ 2 . (76) 
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Using these results in Eq. (l70|) and changing to p a g, one obtains the secular DME in the diagonal basis in the form 

d 



-J t PaP = - ( iuJ a(3 + f p a0 + 8 a p ^ T aa ,p a , a ,. (77) 

In Eq. (|72|) . T aa i is the rate of quantum transitions a' — > a in the small system, accompanied by an appropriate 
transition in the bath so that the total energy is conserved. To the contrary, T a i a in Eq. (|74[) is the rate of quantum 
transitions a — > a' in the small system. One can relate both rates as 

r , - — V p-e*,/v<bT) S , _ , E _ E ,\\v , ,\ 2 

L a' a ' — j-^^ / J c u £-ot f \ ^vd *->v3' ) \ v a>V2,a>'z&' \ 

= V e (- B »'+ e °- £ «')/(fc B T )( j ( _ + ^ _ |y |2 (7g) 

or 

r a , a = e ^-^)/^ T )r Qa ,. (79) 

This is the so-called detailed-balance relation that ensures that asymptotically the small system reaches the thermal 
equilibrium described by Eq. (|67p . If e a < e a > , then at low temperatures the rate of transitions a — > a' (with 
increasing energy) is exponentially small. Note that the transition rates Y a 'a and T aa > correspond to the Fermi golden 
rule. 

The quantity T a p of Eq. (|76[) is the dephasing rate that turns to zero for the diagonal elements of the density matrix 

n a = p aa , (80) 

the populations of states a. The dephasing rate is not related to any transitions of the small system. Its origin is 
modulating its transition frequencies uj a p by fluctuations of the bath. 

One can see that in the decoupled equations for nondiagonal terms of the density matrix a ^ (3 in Eq. ([77)) there 
are only outgoing terms, so that the nondiagonal terms tend to zero asymptotically. The diagonal terms of the DME 
satisfy the system of rate equations 

-^n a = {T aa >n a > - F a r a n a ) . (81) 

The asymptotic solution satisfies 

n , r / e -£ Q '/( fc BT) 

n a ~ T aa , ~ e -"«/<** r > 1 ' 

that corresponds to the thermal equilibrium. The equation for nondiagonal elements can be written in the form 



d_ 

df 

with T a p given by Eq. ([75)) . 



P a p = ~ + f aj gJ p a (83) 
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E. Analysis of the non-secular DME 



We have seen above that within the secular approximation nondiagonal DM elements are decoupled from diagonal 
ones and oscillate with decay independently of each other. In the full DME of Eq. (f63|) . nondiagonal DM elements are 
coupled to the diagonal elements. Nevertheless, they approach zero at equilibrium, in spite of diagonal elements being 
nonzero. This points at an interesting feature of the full DME that has to be worked out in more detail. Separating 
diagonal and nondiagonal elements in the relaxation term of Eq. (|63|) . one obtains 

= -iUaPPap + Rgp,a'0' (l ~ &a'P') P a <0> + ^ Ra(3,a'a' P a > a > 
a'/3' a' 

= -iUaPPu/3 + E RaP,a>f3> (l ~ #a'/3') P a '0' 
a'/3' 

+ E | - E e- E ™ /{kBT) 5 ( e/ j - e 1 + E m - E^)V a ^ 7 ^V^,^p^ 

k vdzu 1 \ 7 

-^T e~ E ™ /{kBT) 5(£ -e +E - E AV ,V , a o 

7 

+ E e~ s »' /(feBT) (S (e -e-y + E^- E w .) V^.^V^.^p^ 

7 

+ e- E ~'' {kBT ^5 [e a -e 7 + E m - E„,) V^^V^^p^ I . (84) 

7 J 

With the use of the energy conservation this can be rewritten as 

^PaP = ~ luJ apP a f3 + ^ R a f3,a'0' (l ~ <$a'/3') P q '£' 
a'/3' 

flZ J J v otvj i i yvJ v •yzu -,pvj 

7 zuzu' 

x ( £/3 - e 7 + £ ro - £ ro ,) (p 77 - e ^-^)/(*BT) Pw j 

+ £ (e Q - e 7 + £ ro - £ ro ,) (p 77 - e (e °- £ - )/(fcsT V Qa )] ■ (85) 

One can see that as the diagonal DM elements approach their equilibrium values, they cease to drive nondiagonal 
elements because of the detailed balance relation. Thus the equilibrium solution of the full DME is Eq. (l67l) . 



F. Semi-secular approximation 



While the secular approximation neglects the interaction between diagonal and slow nondiagonal DM elements, 
the full non-secular formalism involves a big ./V 2 x iV 2 matrix that has to be diagonalized. In important particular 
cases such as thermal activation over a barrier or tunneling, the eigenvalues of the DM span a broad range from very 
fast to very slow, the latter being of a primary importance in relaxation. Because of this, one has to do numerical 
calculations with increased precision that makes them very slow. This difficulty can be overcome with the help of 
the semi-secular approximation that considers coupled equations for diagonal and slow nondiagonal DM elements 
plus decoupled equations for the fast DM elements. The easiest way to implement it is to include the diagonal and 
subdiagonal terms p aa and p a Q±1 into the slow group, because in most situations there are only two levels that 
come close to each other, making p a Q+1 or p a a _ 1 slow. Implementation of the semi-secular DME in the case of 

time-independent H s will be done in Sec. IIII A 31 

G. Transformation to the natural basis 

One can transform the DME, Eq. ([77)) to the natural basis using Eqs. (|2T)|) and (fT^)) in the form 

Pmn = (VUPlVO = ^2{^m\X a )P a p{Xp\^n)- ( 86 ) 

a/3 
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The inverse transformation has the form 

Pa/3 = (XalPlXp) = J2(X a \^ m )Pmn(^n\Xl3) 



(87) 



The general DME in the diagonal basis, Eq. (|63|) , can be transformed to the natural basis with the help of Eq. ([86)) 
as follows 



dt 



Pmn = ^2{^m\X a )-^P a p(Xp\lPn) 



This yields 



a/3 



ml'X.a) \ a - £/9) PapiXpllPn) + ^(^mlXa) R &0,c*0' Pec' P'iX^v 



a/3 



a/3 



a'/3' 



T Y^r,i\X a )£ a {xM m ')P m 'n'{^n'\Xf3){X(j\i'n) 



m'n' a/3 



"% Y^^\X a ){xM m ')P m 'n'{^n'\Xl3)£p{Xp\^n) 
m'n' a/3 

- Y ^WmlXa) ^ R c'P,ct'P'{Xa>\i>m')Pm<n'(' i l J n'\X0>){X0\i>n)- 
m'n' et/3 a' j3' 



^Pmn — y j ^-s, mm' 9 m'n ^ ^ ] Pmn f ^-s,n'n + ^ ] Rmn,m' V P m > n > > 



(89) 



where -ff. 



ft 



Vw ) and 



Rmn.m'n' — ^ ^ ^mn.m'n' ;a(3,<x f i Rql(3,q.' (3 1 3 



where 



mn.m'n 



;a/3,a'/3' = (V'm I Xa) (X/3 I V'n) (Xa' I V>m' ) (VV I X/3' 



(90) 



(91) 



Additionally, the matrix elements in the Fermi-golden-rule transition rate, Eq. (|72p . that are defined with respect 
to the diagonal basis, can be expressed through those with respect to the natural basis. Similarly to Eq. (|87|) one 
obtains 



Varo.a'ra' — ^ (Xal'0m)^7»'gP,m'n; / {^m' IXa')- 



(92) 
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III. TIME-DEPENDENT PROBLEMS 



In this section wc consider the DME and its solution in three important cases: i) free evolution for time-independent 
H s ; ii) fast resonance perturbation and iii) periodic or nonperiodic slow perturbation. In the second case it is sufficient 
to keep the perturbation in the conservative part of the DME only, while in the third case modification of the relaxation 
terms is required, as well. In all these cases the solution can be obtained by matrix algebra. To the contrast, problems 
with a large temporal change of H s cannot be solved by matrix algebra. The secular, non-secular, and semi-secular 
versions of the DME yield the same results except for the case of anomalously close energy levels (e.g., tunnel split 
levels). In the latter case the semi-secular DME is preferred, while in general the fastest and easiest secular DME is 
the best choice. 



A. Free evolution 



1. Non-secular DME 



The solution of time- independent DME, Eq. (|63|) . that can be rewritten as 

= ®af3,a'/3'P a <l3'y $a(3,a>l3' = —lUapSaa 1 & ftp + R a f3 ,a' /3' , (93) 

a'/3' 

is a linear combination of time exponentials with exponents being eigenvalues of the matrix <I> building the DME. To 
bring the DME into a standard form, it is convenient to introduce the compound index a defined by 

a = a + N(P-l), (94) 

where N is the size of the density matrix, i.e., a, /3 = 1, 2, . . . , N. Then a = 1, 2, ... , N 2 . Inversion of Eq. yields 

a = 1 + ATrac ( , = 1 + Int ( . (95) 



N J \ N 

With the index a, the density matrix p a p becomes a vector with the components p a while $ a ^ a /^ becomes a matrix 
with the elements <I> aa ' : 

a' 

The eigenvalue problem for the DME can be written as 

* ■ R /1 = -A^R^, A„ = r p + m^, (97) 

where is the right eigenvector corresponding to the eigenvalue A^ and fj, = 1,2, ...,(25 + l) 2 . Since $ is a 
non-Hermitean matrix, right eigenvectors differ from left eigenvectors that satisfy L^-fl? = — A^L^. Left and right 
eigenvectors satisfy orthonormality and completeness relations 

(98) 

In general, L^ and R^ are not Hermitean conjugate, see, e.g., Eq. (|103|) . All real parts of the eigenvalues are positive, 
> 0. There are N purely real eigenvalues, one of which is zero and corresponds to thermal equilibrium. We assign 
the zero eigenvalue the index /x = 1. Complex eigenvalues occur in complex conjugate pairs. 
The solution of the DME with the initial conditions can be written in the form 

p(i) = ^R / ,e- A - t L /J -p(0). (99) 

The fully vectorized form of this equation is 

p(t) =E- W^-E" 1 -p(0), (100) 
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where E is right-eigenvector matrix composed of all eigenvectors R M standing vertically, E _1 is the left-eigenvector 
matrix, composed of all left eigenvectors lying horizontally, and W(£) is the diagonal matrix with the elements e~ Afl *. 
In fact, 

E • W(t)-E _1 = exp . (101) 

The asymptotic value of pit) is described by the zero eigenvalue Ai = 0, 

p(oc) = Ri (Li • p(0)) . (102) 

Here p(oo) = p cq should be satisfied, where p eq follows from Eq. (|67[). and this result should be independent of p(0). 
Thus one concludes that 

Lla = 5 a (a)l3(a), Rla = ]T CX P y'feT^ ^«(<0/3(<0 ' ( 103 ) 

where a (a) and /3(a) are given by Eq. (f95|) . This means that Li is related to the normalization of the DM while Ri 
contains the information about the equilibrium state. One obtains 

Li • p(0) = £W o (0) =^5^(0) = ]Tp QQ (0) = 1 (104) 

a a/3 a 



and p(oo) = Ri — p cq , as it should be. Note that Ri and Li satisfy the orthonormality condition in Eq. 

Y^LiaRia = T^^cxp f-^j S a{a)l3{a) = — ^ exp (-7^) <W = 1. (105) 



Z s <^ r V k B Tj awPW Z s ^ ' \ k B T, 

a a x * a/3 x 

The time dependence of any physical quantity A is given by Eq. Q that can be rewritten in the form 

A(t) = ]T A p a p aP (t) = J2 AaPjt), (106) 

a/3 a 

where 

A Q = Ap [a)a(a)l A 0a = (p Act). (107) 
Writing Eq. I|106|) in the vector form as A(t) = A • pit) one obtains 

A(t) =^A-R Al e- A - t L, I -p(0) (108) 

A* 

or, in the fully vectorized form, 

A(t) = A E W(t)-E _1 -p(0). (109) 

Since the time dependence of observables in the course of evolution of the density matrix is described by more than 
one exponential, one needs an appropriate definition of the relaxation rate or relaxation time. A convenient way is to 
use the integral relaxation time defined as the area under the relaxation curve 

_ j~dt [#)-iH 

int " A(0)-A(oo) • 1 Uj 

One can check that in the case of a single exponential, A(t) = A(oo) + [A(0) — A(oo)] e~ rt , the result is r- m t = 1/r. 
From Eq. (|108[) one obtains 

A(t) - A(oo) = ]T A ■ R ^ e " A "' L M ' P(0) (HI) 



and thus 



E^ 2 (a-R^a^(Lm-p(o)) 

T *nt = „at2 ; . — - ; ~ — — ■ (ll 2 ) 



Jfj,=2 "/V iy /i V-^M 

This formula cannot be fully vectorized since summation skips the static eigenvalue p = 1 
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2. Secular DME 

Within the secular approximation, one has to consider the dynamics of diagonal and nondiagonal components of 
the density matrix separately. The former is described by Eqs. (|96j) - (|102j) where the vector p is replaced by the vector 
of the diagonal components n = {n a } = {p aa } and the (N) x (N) matrix 3? is replaced by the (N) x (N) matrix 
<1> SCC having matrix elements 

C- = (i - <W)r Qa ' - <W ^r 7Q , (113) 

7 

as follows from Eqs. (|TT[) or (|8Tj) . All eigenvalues of <& scc are positive reals, except for one zero eigenvalue, Ai = 0. 
Eq. (|103|) becomes simply 

L la = 1, R la = i- cxp (—j^f) ■ ( 114 ) 

If the initial condition is a diagonal matrix, the non-diagonal elements do not arise dynamically and hence they can 
be dropped. Then the time dependence of any quantity A is described by 

A(t) = A aa n a {i) = A E W^-E" 1 • n(0), (115) 

Oi 

c.f. Eqs. (|106[) - (|109p . For the integral relaxation time in the case of a purely diagonal evolution one obtains 

S^(A-R,)A,'(L,-n(0)) 

Tint = t^-Tr , (116 

EjU(A-R M )(L M .n(0)) 

c.f. Eq. fTT2]) . 

If the initial state is a non-diagonal density matrix, one has to add the corresponding trivial terms following from 
Eq. HMD, 

A(t) = A E W^-E- 1 • n(0) + ^ A^ Paf3 (0)e~ (l ^ + ^' 3)t . (117) 

Then Eq. (| 1 16[) is generalized to 

EjLa (A • R M ) A; 1 (L„ • n(0)) + £ a ^ A 0a p a e(O)(iu af} + f^)" 1 

Tint = m • (H°j 

E^ =2 (A • R„) (L„ • n(0)) + A PaPap (0) 

Obviously this expression is real. 

3. Semi-secular DME 

Within the semisecular approximation introduced in Sec. Ill Ft the slow group being formed by diagonal and 
subdiagonal DM elements, \a — j3\ < 1, the equations of motion for the latter have the form 

~M P <*& = X! ( $ «/3;a'.«'-lPa',a'-l + ®af3; a > a' P a > a > + $af}; a > ,a>+lP a > , a '+l) ( 119 ) 



that is a subset of Eq. ([93]) . In labeling matrix elements, one can introduce the compound index 

a=2(a-l) + /3. (120) 
Here a = 1, 2, . . . , TV and /3 = a — 1, a, a + 1, so that a takes the values a — 1, ... , 3N — 2. In terms of a one has 

a = l + Int(!), /3 = a - 1 + 3Frac (J) , (121) 
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and Eq. (|119[) can be rewritten as 



!/?f w =EW^° W < (122) 



where $ aa / = $ a (a),p(a)jatar),p(a')- 

The solution of Eq. (|122p is similar to that of Eq. On the other hand, there are uncoupled DM elemens with 

| a — /3\ > 2, like in the secular approximation. Instead of Eq. I|117|) one has 

A(t) = A-E-W(t)-E- 1 -p slow (0)+ ^«^(°) e_(lU ° ?+f ° 3)i (!23) 

Q-/3|>1 



and instead of Eq. (] 1 1 8[) one has 



E^ 2 " 2 (A • R M ) A; 1 (L„ • p slow (0)) + E\ a -m>i Ap aPap {0)^ afj + f^)" 1 



(124) 



(A ■ R M ) (I> ■ p^-(0)) + E\ a -p\>x Ap a p aP (0) 
B. Resonant perturbation 

Let H s contain a periodic perturbation 

V(t) = V Q e- iut + Vje iut (125) 

with the frequency lu > close to one of transition frequences uj m i > 0. For the resonance to occur, the latter should 
exceed the relaxation rates, so that one can use the secular DME, Eq. (177]) . Since Vq is small, one can disregard it in 
the relaxation terms and use the diagonal basis corresponding to the time- independent part of H s . Combining Eqs. 
(|52| and ([77]) , one obtains 



0' 



(iuJafi + T a p) p a(i + 6 a p 

E r ^P 77 - (126) 



Here one should drop all nondiagonal DM elements but p„„, and p^,^ since the former are not excited by the resonance 
perturbation and vanish with time. As a result, one obtains a coupled system of equations for the diagonal elements 
Paa an d the elements p m i and p^^ = p*^ of the form 



d i 

J t P m = ~\ {e-^Vo^p^ - e^p w (Vj) ^ + £ (r, 7 p 17 - T lvP J 



j t Pn'n' = -\ ( elUt *V ^ e-^p^K,,™') + ]T (r,, 7 p 77 - T^p^,) , (127) 

plus Eq. (|8T|) for all diagonal elements with a ^ J], rj '. Also terms such as e lut (Vq) p v / v ~ e l ( w+a; '»''')* in the second 

equation have been dropped. Neglecting such terms that oscillate out is the so-called rotating-wave approximation 
that is similar to the secular approximation. It is convenient to introduce the slow nondiagonal elements p m , and p^^ 
via 

Priri' = Pr)7]' e ' Prj'r] ~ Pr]'r) e ■ (128) 
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With the use of Vq 1 = V* , the DME acquires the form 



'/'»/ 



d . 
d 

dt P 
d 

dt P * 
d 

dt P ' 



Im (p nri > V Q riri ,) + (r^ 7 p 77 — T^p^ 



Pr/rj' 



^ I m (ftp)' K),?)??' ) + J]] (pV^Pfl ^~fV' Pri'v' / 

77V 



7#" 



(129) 



Here f,,,,/ is given by Eq. (|75|) . 

In the strong-dephasing case r w ' 3> r 77? /, Vo. the nondiagonal element p r , quickly reaches its quasistationary 
value that can be obtained from the first equation of Eq. (|129[) by setting dp m ,/dt = 0. This results in 



1 Vb,i?i7' (Pgv - Prrn) 

Prtr}' t ~ • 

n w — Lu rm , + iT m , 

Substituting this into the equations for the diagonal DM elements, one obtains 

2 



(130) 



dt Prm 



tit 9 "'* 



T2 \Pr)'n' ~ Pnn) 7 \2 . f, 2 ^ (JVr^VY _ ^lvPr,v) 



0,777;' 



/; 2 



\Prj'r]' Prjrj) 



2F / 



(r i; / 7 p 77 ^■yrj'P v ' 7 y) 



d_ 

dt' 



^ ] (r a7 p 77 -T 7Q p QQ ) , a i],r] . 



(131) 



Note that the absorbed power of the field acting on the system is given by 

2 



P = flLOrjr,' 



Vo. 



</'/' 



ti 2 



{^Prj'rj' Prjrj) 



2F / 



(uj - uj w ) 2 + f 2 m , ' 



(132) 



since every transition \rf) — > |?7) absorbs the quant of energy equal to hw m * . Eq. (|13ip can be simplified by introducing 
the experimentally measured power P and thus eliminating the matrix element Vo, W ' and the dephasing rate T vv r. 
Still the equations are rather complicated so that one cannot find a simple general solution even for the stationary 
state. 

Let us consider at first the model with only two levels, 77 and rj'. Here with p m = 1 — p and the induced-transition 
rate 



A = 





2 




2F / 




(W - Urpf) 2 + f ^ 



(133) 



Eq. ()131|1 reduces to a single equation 



d_ 

Jt 1 



A \Ptfiy Pr]rj) + ^V'V Prjrj ^VV' Pr/'ri' 

= - A (2P V > V > T v'n ~ (Tri'l + T vv')Pv'v' 

= A + T r y rl - (2 A + T v > n + T m i) p, q , ri , 

that describes relaxation with the rate 2A + r n > v + r ?))) < towards the stationary state described by 

A + T,y rl 



Pr/'rj' 



2A -|- r r ^^ -t- 



^rjrj' ^rj'rj 6Xp 



(134) 



(135) 
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for the ground-state population. The limiting cases of this formula are 



Prj'rj' 



„(eq) 
rj'rj' 



1 + exp 
1/2, 



, A = 
A — > oo, 



(136) 



as it should be. For the population difference from Eq. (|135p one obtains 

Prj'rj' Pr\r\ ^Pti'tj' ^ 



r / — r / 

*7 7? r l r l 



(eq) _ (eq) 



2 a + ry„ + v m , i + 2A/(rv„ + r w ) ' 



(137) 



One can see that the population difference is reduced by the resonance perturbation. The absorbed power is then 
given by 



fllOr 



P = fiw w (p v , ri , ~ PrjV ) A = — 



2A/(T^+r w ) ' 



(138) 



It increases with A and reaches an asymptotic maximal value. 

To relate the populations directly to the measured absorbed power P, it is most convenient to rewrite Eq. (|131[) for 
the two-level system in the form 



d_ 



This equation has a stationary solution 



P 



P 



J rT 1 y v p rm T nn i p^i^i 

+ ^ rj'rj — rj'rj + Tjjjj/) p^i^i ■ 



_ T^ v - P/(hlO m >) _ ( eq ) P/(huJ m ,) 
Prj'rj' i i , P' 



1 rj'rj i - 1 - 7777' 



(139) 



(140) 



I77' rj ~t~ r???/ 

Note that the maximal absorbed power corresponds to the saturation, p , , — p „, i.e., p rjlri , = 1/2, wherefrom follows 



P, 



o V Tj'jj ^ rjrj' J n 1 



v n 



1 — exp 



tv-^rjij' 



(141) 



that also can be obtained from Eq. (|138[) . 

Consider now a physical quantity A that has the expectation value given by Eq. ([4]). The contributions from the 
nondiagonal DM elements oscillate in time and thus average out. The result is due to the diagonal elements only, 



A — A r j r jP vv + Ajjijf p„,„i — Ay m + (Ajjijy A rm ) Pjyjji 

With the help of Eq. (|140|) it can be rewritten as 

A = A^ + AA, 

where 

A^rj' + A m exp ^ k^r^j 



^( c q) _ _|_ {Arfrf A m ) p^, — 



1 + expl-^ 



is the equilibrium value and the deviation from the equilibrium is given by 

A , , — A P 

s^-rj'rj' s%7jrj 1 



A A = -■ 



^r/'rj ^rjrj' fyjJrjrj 1 

Thus the total relaxation rate between the states \rf) and \rf) can be expressed by the formula 

A — A , , P 

finn fin' n' 1 



r = r / + r > = 

rj'rj ~ rjrj 



^vv__2zlJi_ 

AA ftw m ' 



(142) 



(143) 



(144) 



(145) 



(146) 
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through the measured P and A A. 

If the system has more than two levels but the levels rj and rj' are the lowest two levels and the temperature is 
low, so that the populations of all other levels are small, one can still consider an effective two-level model in which 
relaxation between r\ and 77' can be assisted by the upper levels. Instead of the direct relaxation \rj) — ► \rj') (that can 
have a small rate) the system can be thermally excited from \rj) to some high level \a) and then fall down to \r]') . 
This is the Orbach mechanism that has the characteristic Arrhenius temperature dependence of the rate. Using the 
effective two-level model is justified by the fact that the upper levels do not contribute to physical quantities A at low 
temperatures. 



C. Linear response 

In this section we consider a small harmonic perturbation of Eq. (|125|) acting on the small system, similarly to the 
preceding section. However, the frequency to of the perturbation does not need to be close to the resonance with any 
transition w a p. If to <C uj a /3 for all a, /3, the response of the small system is due to the relaxation and it depends on 
the relation between lu and the relaxation rate T. To make an account of this effect, one has to include V(t) into the 
relaxation terms of the density-matrix equation. Temporal change of V(t) changes the instantaneous equilibrium to 
which the system relaxes that gives rise to a dissipative dynamics. On the other hand, one has to include V(t) into 
the conservative term of the DME as well, where it can work as a resonance perturbation. The secular approximation 
will be used below for simplicity. 

Using the adiabatic basis of Eq. (|56)) and combining Eq. (|58|) with Eqs. (|8Tj) and ([83)1 , one can write the secular 
DME in the form 

j t P aa = Y l*r) p t° + Pcnbc*, l*a» + Y Ta ' a (e {£a '~ ea)/{kBT) P a 'a> ~ P aa ) (147) 

7 a'^a 

for diagonal terms and 

^Pais = Y (ix a |x 7 ) P 1 p + P ai {x 1 \xp)) - (iu a /3 + fa/3) P a p (148) 

7 

for nondiagonal terms. If time derivatives of the adiabatic states are small, the density matrix in the adiabatic basis 
p a p is close to its instantaneous quasiequilibrium form given by Eq. (|67p . The time dependence of the adiabatic states 
drives the DM out of the instantaneous equilibrim, p a/3 lagging behind time-dependent p^p. On the other hand, since 

V(t) is a small perturbation, p^ deviates from the full equilibrium p„ C g eq in the absense of Vit) at linear order in 
V{t). Thus the deviation from the full equilibrium, 

SPa^Pap-C^ (149) 

is also linear in V(t). Since the difference e ty£a '^ £a ^^ kBT ^ p a , a i — p aa and the nondiagonal elements p a/3 are small, one 
does not have to expand uj a and the relaxation rates. The only term to expand is e^ £ °''^ £ ''^^ kuT \ where one can use, 
at linear order in V(t), 

£a {t) = £ (°) + 5e a {t) = \V{t)\ x (°>) = V(t)M. (150) 



Thus 



e (e a ,-e«)/(k B T) s e ( E l?- e ««>)/(t B T) L + V^L' ~ V(t)^ \ ^ 



The driving term expands as follows: 

(&« Wi) PiU + Po n {x 1 \xp)) 

7 

= Y{^\^)p ( T+p { «V(x,\xe)) 



(X a \X P ) Pfr + P ( ^(x a \X P ) = (X a \xp) (Pi°2 cq - pfP) , (152) 
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where (x a \Xp) = &a(3 was used. One can see that in the diagonal equations this term disappears. Thus for the 
diagonal equations the linearized DME has the form 



dt 



.( e W- e g))/ffa T ) .(0)eq VWSL ~ ^(*)«» 
"' "a' a' rp 



(0) 



a'^a 



or 



where 



J act — y / 1 a' a 



at' 

Linearization of nondiagonal equations yields 

= - ( w a(3 + f a/j) Pa/3 + /a/3, 

where 



(153) 



(154) 



(155) 



(156) 



(157) 



One can see that within the secular approximation the driving term is conservative in the nondiagonal equations and 
dissipative in the diagonal equations. 

The density matrix p a g in the equations above is still defined with respect to the adiabatic basis |x Q ) , although the 

relaxation term is already expanded around the unperturbed states 



Xa^ ) • Let us now make transformation to p®l 



defined with respect to the unperturbed diagonal basis 
the form 

,(°) 



) • The relations between the two density matrices have 

(158) 



Pal = E^ \Xa>) Pa'/3'(Xf3 
a' 13' 



(0) 



c.f. Eq. IJSBJ), and 



Pa/3 = E < X * X »') Pa'P'^P Up) 
a'P' 



(159) 



Differentiating the first equation over time one obtains, at linear order in V(t), 



dt' 



a'p' 



Pal = ^(X^lXa^Pa'P'iXp'lxp j + ^iX^lXa^Pa'p'iXp' Xp ) + ^{X ( a \X a ') K' f3' {Xp 



a' 13' 



a'f3' 



(0) 

X(3 



E<xi 0) \x a ')& q +T,pT^' \xf) + E<^ 0) \x a ') Pa ',{x P ' \xf 



0' 



a'/3' 



-<xi 0) \x ) (pi°2 cq - pJT) + E<xi 0) \x«) fork? 



a' S3' 



(0) 

Xf3 



(160) 



For a = the first term in this equation disappears and for a ^ f3 it cancels a similar term in Eq. (|156p . Thus Eq. 
(|156| with a j3 transforms as 



(o) 

X/3 



-J t P a j3 = - E^ \X-oi) (iUa'l3> + f a'/?') Pa'P'iX/3' 
a'f3' 

= - E E ^" 0) l^a') ( iw «'/3' +fa'/3') (Xa' Xa" ) P«/ '/3» <X^» |X/3') <X/3' 



a'/3' a" /3" 



(0) 

x^ 
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In this expression the projectors such as (\a^ \x a ') are linear in V(t), if the indices do not coincide. Thus there are 
two types of contributions: (i) All indices pairwise coincide and thus pf) l/3 „ — * p^p ~ ^(*) or (ii) a" = f3" so that 
Pa"/3" ~ * ^a°''a" an( ^ onc °^ tnc P a i r °^ indices in projectors do not coincide. This yields 



- O (0) 
dt & 



,(°) 



=,(0 



J a'ff ~ x a 



)) \ (0) 
J ^a/3 



£<xi 0) lxa<> (few + *W) p^<x2 5 |X^> (X,' 



a'/3' 



E(xi 0) IXa') (*^a'/3' +fa'/3') (Xa' 



a'/3' 



v (0)\ „(0)eq, 
X/3' //>/3'/3' \X/3' 



(0) 
*/3 



(0) 

4 



(°) ^ J°) 



-<xi 0) \x P ) tfvP^ - r aaP ^( Xa |x? } ) - + f afJ ) P ^(x^ \x ) 



(0) 



+ <Xa|x? ) )- (0,< " 1 



(161) 



i.e., 



A JO) 

dt 9 ^ 



-(lJ a) +f (0) W 0) 

-2( X £°) |x,) f w pS cq - 2 <Xa |x? } ) f aaP ^ 

-(toa/J + f^) ((xi° ) |x /3 )pi°2 Cq + (X Q 



The small projections (xi°' ) |X/j) can b e calculated with the help of the perturbative expansion 



IXa) 



^ y Xa' 



(0)\ 



,(0) 



JO) JO) 
fc a — 



(162) 



(163) 



that yields 



(xi 0) |x„> 



V(t) 



(0) 



JO) _ JO) 

& ta 



(0) 

a/3 



(Xo 



(0) 



Y (0) 
Xa 



K(t) 



(0) 

X/8 



.(0) JO) 



(0) 
a/3 



Thus the DME for nondiagonal components has the form 

- o (0) --^ (0) + f (0) W 0) + f 



(0) 

a/3 ' 



where 



JO) _ 

/a/3 — 



^(0)' 



1 + 



(0)cq _ (0)eq 
aa P/3/3 



') 



f, (0)eq f, 
r a aPaa ~ T^p 



(0)cq 
5/3 



,(0) 



(0) 

a/3 



J a/3/ ^a/3 

As the secular approximation implies T a p <C u) a p, relaxation terms in this equation can be neglected, 

,(0)eq _ „(o)oq^ 



(164) 



(165) 



(166) 



(167) 



The conservative part of this driving term could actually be transformed exactly using the relation between the two 
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bases, Eqs. (l52l and (|58|) . Transformation of the diagonal Eq. (|154[) is as follows: 
d 



,(o) 



dt' 



y(0) 



a" 



E<* 

a IXa 



(Xa' 



/(°) 



E ( r -' E x^)pi°V(4" l^)-r^ Q E (x. xS^p^VW IxJ ] +/£ 

a"/3" a"/3" 



V (r ,o {0) -r - o (0)N l + f 



(0) 
a a ■ 



Since there are no corrections to r aQ / linear in V(i), this finally yields 



f(0) 



where 



(0) 



V r (0 ,> 

/ j a' a 



(o) V{tt) a ,-V{ttl 



(168) 



(169) 



(170) 



One can see that this transformation was trivial. In the sequel we will drop the superscript (0) . 

Let us now solve Eqs. (|165p - (|170| . Similarly to Sec. IIII A2( one can write Eq. (|169|) in the vectorized form 



dt 



-Sn = * S0C Jn + f diag , 



(171) 



where (<5n) Q = Sp aa and the elements of «fr scc are given by Eq. (|113|) . Since f contains positively- and negatively- 
rotating terms, the stationary solution of this equation has the form 



± iu6n {±) = $ sec -<5n (±) + f dia s>(±) 



with 



e -e a /(k B T) 

f(±) — \ p , y a'a 

J act g / j 1 a' a 



V (±) _ y(±) 



(172) 



(173) 



Here V$*\ = Vo,« 7 and 



v 0,a~f — \ V 

\ /a. 

eigenvectors defined by an equation similar to Eq. (|97p . 

^ (±) = E^ ±)R 



V$„ a . The solution of this equation can be expanded over the set of right 



(174) 



Inserting this into Eq. ()172|1 . multiplying from left by the left eigenvector L„ and using orthogonality in Eq. |98|) . 
one obtains 



and 



± iwCW = -A v Ci^ + L„ ■ f diag ' (±) 
Now the final result for the populations is 



A,, ± ILU 



Afj, ± iuj 



(175) 



(176) 



(177) 
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At equilibrium, u — 0, this expression should reduce to the static result 

T . f diag,(±) e -e a /(k B T) / g §Zt 



6 



n, 



E E ^ = ( + f J ^ 178 ) 



that can be proven to satisfy Eq. (|169p in the static case. Using Eq. (|173[) and ^ ± V^q" 1 = Se a , one obtains the 
identity 

"X7^ a " a ' k^f MQ '""^r~X 1 j 

that should be satisfied by the matrix solution and can be used for checking. Nondiagonal components of the DME 
satisfy the equations 

± iwSpW = - (iu> a p + f afj ) 5p ( ^ + f { a f , (180) 

where 

These equations have the solution 

f (±) 

fip$ = * ; a0 . . ■ (182) 

For a physical quantity A the linear response has the form 

A(t) = e luJt A^ + \oj) + e-^A^iw), (183) 

where 



A(± = » (A.R,)(L,.^W) + y A fia f% 



(184) 



c.f. Ecj. pi8[) . Here (A) = A aa . The zero eigenvalue, /i = 1 and Aj = 0, does not make a contribution to this 
formula since Li given by Eq. (|1 14[) is orthogonal to f dia S'(±) defined by Eq. (|173[) . Physically relevant are real and 
imaginary parts of the linear response 

A'(u) = ReA^(iv) +ReA(-\uj) 

A"(u) = -Im#)H+Imi(-'(w). (185) 
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IV. APPLICATION TO MOLECULAR MAGNETS 
A. The model 

Let us consider, as an example, a magnetic molecule (MM) described as a large spin S with the Hamiltonian 

H S =H S = H A + Hz, (186) 

where Ha is the crystal-field Hamitonian and Hz is the Zccman Hamiltonian, Hz — —gfJ, B H ■ S. The crystal field Ha 
is usually dominated by the uniaxial anisotropy, 

H A = -DS 2 Z - BS^ + H' A , (187) 

whereas H' A contains smaller terms that do not commute with S z . In the absence of the latter and of the transverse 
field, the eigenstates of the spin are |m) , m = —S,...,S. The energy levels of the magnetic molecule are given by 

e m = —Dm 2 — Bm A — grriBH z m. (188) 

The transition frequency for a pair of levels is 

hw mm i = e m — £ m ' = — (m — m!) { (m + to') \D + (m 2 + to' 2 ) B] + gmsH^ (189) 

Condition Hu) mm i — for to ^ to' (levels on different sides of the barrier created by the uniaxial anisotropy) defines 
the resonance values of the longitudinal field H z . For the generic model of molecular magnets with B = the latter 
are given by 

gm B H z =kD, k = 0, ±1, ±2, . . . (190) 
For these fields all levels in the right well 

m = —m — k (191) 

are at resonance with the corresponding levels in the left well m < 0. For the realistic model with B > 0, the field 
creating resonances between low-lying levels with large m 2 + ml 2 is greater than the resonance fields for high levels. 
Transverse anisotropy and transverse field H x that enter H' A result in the tunneling under the barrier and tunneling 
splitting of the resonant levels m, to'. 

For Mn i2 used in illustrations below we adopt the values D/Ub — 0.548 K and B/ks = 1.1 X 10~ 3 K (Refs. @,[1,[B|) 
that makes up the barrier of 66 K. H' A in Eq. (|187j) can contain second- and fourth order transverse anisotropy, 

H' A = E(Sl- S 2 ) +C{Sl + St). (192) 

For Mn 12 C/k B = 3x 10~ 5 K (Refs. |,||), whereas E = in the ideal case because of the tetragonal symmetry of 
the crystal. However, it was shown [fjj that local molecular environments of Mni2 molecules have a two-fold symmetry 
and rotated by 90° for different molecules. Although on average the four-fold symmetry of the crystal is preserved, it 
gives rise to nonzero E that will be set to E/k-Q — 2.5 x 10~ 3 K. 

The spin Hamiltonian of molecular magnets can be easily numerically diagonalized to yield eigenstates \a) and 
transition frequences u/ a p, a,/3= l,..,2S + 1. The terms non-commuting with S z such as H' A and the transverse 
magnetic field cause hybridization of the spin states in the two wells that leads to spin tunneling. 

B. Spin-phonon interaction 

The magnetic molecule is embedded in the elastic matrix described by the harmonic-phonon Hamiltonian 

H h = H ph = ^ ^kAak A a k A- (193) 

kA 

Describing the spin-phonon interaction, we will follow the approach developed in Refs. @, @] that allows to avoid using 
unknown spin-phonon coupling constants and to greatly simplify the formalism. As a magnetic molecule is more rigid 
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than its ligand environment, a good approximation is to consider this molecule rotated by transverse phonons without 
distortion of its crystal field. This leads to the spin-phonon interaction 



V = ff s -p h = RHaR- 1 -H a , R = e-*-** 
where 5<fi is a small rotation angle given by 



1 



V x u(r), 



(194) 



(195) 



u(r) being the lattice displacement due to phonons. Expanding Eq. (| 1 94[) up to the second order in Scf> components 
yields 



V = V {1) + V (2) 



where 



and 



Ha, S 



6<f> 



y( 2 ) = — 
2! 



Ha, <?£ 



, St 



x,y,z 



with summation over repeated indices. We use canonical quantization of phonons, 

„zk-r 



U 



2MN ^ 

kA 



(a k A + a T _ kA ) 



(196) 
(197) 

(198) 
(199) 



where M is the mass of the unit cell, N is the number of cells in the crystal, ekA are unit polarization vectors, 
A = t\,t2,l denotes polarization, and u>k\ = v\k is the phonon frequency. The operator 5<fi that follows from Eq. 
(|195[) is given by 



5d> = - 
v 2 



h v^[ikx e kA ] e 4k r 
2MiV^ ,akA 

KA v 



*-kA 



(200) 



Only transverse phonons, etA-Lk, survive in this formula. Whereas is linear in phonon operators and describes 
direct phonon processes, is quadratic and describes Raman processes. Relaxation rates due to Raman processes 
are generally much smaller than that due to the direct processes since they are the next order in the spin-phonon 
interaction. However, the rates of direct processes can be small for special reasons, then Raman processes become 
important. Processes of orders higher than Raman can be always neglected. 

It is important that the spin-phonon interaction above does not include any poorly known spin-lattice coupling 
coefficients and it is entirely represented by the crystal field Ha- Moreover, the relaxation terms in the DME can 
be represented in the form that does not explicitly contain Ha, the information about it being absorbed in the spin 
eigenstates |a) and transition frequencies uj a p that can be found by numerical diagonalization of H$. This can be 
achieved either by changing from the laboratory frame to the local lattice frame in which Ha remains constant but 
an effective rotation-generated magnetic field arises 0, H, Q or by manipulating matrix elements of the spin-phonon 



interaction with respect to spin states, 
particular, for V^ 1 ) one can use 



V 



P), Ref. 



Both methods are mathematically equivalent [8]. In 



H A ,S = 



H s - Hz , S = 



H S ,S 



iSxgn B H 



and the fact that |a) are eigenstates of H$ to obtain the spin matrix element 



For one writes [T(| 



-afi 



Ha,S^ ,S^> 



H A ,S 



(3) = ihw aP (a |S| 0) - (a |S| p) x. g/ i B H. 



(201) 



(202) 



Hs — Hz J , S$ 



Hs , Si 



g^B 



H e S 6 



(203) 
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Here the first term can be transformed as follows 



a 



Hs, S? 



= [( £ « - *r) <« 1^1 7) <7 \S e \ p) + (ep - e 7 ) (a |%| 7 ) (7 I S« I /?>] 
7 

=► ^(e a + e/3 -2e 7 )(a|5 4 | 7 )( 7 |%|/3), (204) 



where on the last step the symmetry properties of Eq. (|198|) were used. Finally one obtains 



= i 2 < a 



iT A , S f l , %] p) = -Y,(s a + s f3 - 2e 7 ) (a |^| 7) (7 |% | P) 



(205) 



Eqs. (I202j) and (|205|) provide a great simplification of the formalism, since otherwise one would have to derive different 
forms of the relaxation part of the DME for each particular Ha- 



C. DME for molecular magnets 



1. Secular vs non-secular 



Most of numerical work on molecular magnets used the secular form of the DME, in fact reduced to the system of 
rate equations for the populations. However, tunneling resonances can make the secular approximation invalid. Indeed, 
if two levels a and a' of the small system have very close energies, the density matrix element p aa , is oscillating with 
a very small frequency Lu aa > in the absence of the coupling to the bath, see Eq. (|63p . If Lu aa ' is smaller than the 
relaxation rate between the neigboring energy levels, p aa , does not decouple from the diagonal elements p aa and p a i a /, 
and the secular approximation breaks down. It is easy to demonstrate that the failure of the secular approximation 
at resonance may lead to unphysically high escape rates out of the metastable state. Consider a MM exactly at kth 
resonance with k > and a very small H' A and the transverse field, so that the metastable ground state \—S) is at 
resonance with the excited state \S — k) in the right well. The latter can decay into the lower- lying state \S — k + 1) 
with the rate Ts-k+i,S-k- Since the exact eigenstates at the tunneling resonance are |±) that are linear combinations 
of |— S) and \S — k) (see Sec. IIVE[) . both of these eigenstates are damped with the rate of order Fs-k+i.S-k ( m 

fact, half of it). The secular DME uses rate equations for p ++ , p , etc., and the initial condition spin in the state 

|— S) gives rise to the initial conditions p ++ (0) = p__(0) = 1/2. Since both p , , and p__ relax with a rate of order 
Ts-fc+i.s-fc, the spin quickly leaves the metastable state, even in the case of a vanishing tunnel splitting, A — > 0. 
Indeed, such an unphysical behavior follows from the analytical and numerical solution of the secular DME at weak 

tunneling resonances. In contrast, coupling of p ++ and p to the slow non-diagonal DM elements p^ and p |_ in 

the non-secular DME leads to the physically expected vanishing of the escape rate in the limit A — > and T — > 0. 

Below the non-secular DME, Eq. will be used in the development of the formalism. The secular and semi- 

secular reductions of it can be obtained later. The relaxation tensor R a p. a if3> is a sum of two contributions, 

R a p,a>l3' = R ala'l3' + ^o/L'/J'' ( 206 ) 

that are due to the first- and second-order phonon processes. These contributions will be calculated separately below. 



2. Initial condition for free relaxation 



Let us consider the question of the initial state of the spin in the case of free evolution. In resonance experiments it 
is, typically, the first excited state. Although, practically, in these experiments only a small portion of the population 
is being transferred from the ground state to the excited state, one can consider the system prepared fully in the 
excited state because of the linearity of the DME. Preparing the spin in the metastable energy minimum, one can 
study its thermal activation over the barrier and tunneling under the barrier. In general, it is not easy to find the 
quantum-mechanical state realizing or approximating this classical state, and in the case of a tunneling resonance such 
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a state does not exist. A good practical way to create such an initial condition is to prepare the spin in the coherent 
state |n(#, (p)) pointing in the direction of the metastable minimum found classically. The spin coherent state is given 

by 



\n(d,<p)) = C m \m), 

m=-S 



where 



C m — 
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S- 



1/2 



S+m 



cos-. 



sm ■ 



S—m 



(207) 



(208) 



3. Direct processes 



To compute matrix elements , etc., in Eq. ([51)1 with respect to the phonon bath, one can label the state 



by the numbers of phonons ^ kA = 0, 1, 2, ... in each phonon mode kA 

l</0 = I ...jfkA,...) =>■ WkX)- 



(209) 



In the direct processes, the state \4> m i) is not independent and it differs from |</> ro ) by creation or annihilation of one 
phonon, according to Eq. (|200[) . We will make use of the phonon matrix elements 



M±(k) = (i/ kA ± 1 |<50Ka 
and their conjugates. From Eq. (I200j) one obtains 

1 / 1T~ [ik x e kA ] e ikr 



(210) 



M_(k) = ^ kA - I 
M*(k) = ( ^ kA 



2 V 2MN 



-a kA 



^ kA 



1 / h [ik x e kA ] e ikr 



2 V 2MN 



1 h [-ik x e kA ] e lkr + 
a 



M+(k) 



2 V 2MiV 
^ kA + 1 



1 h [—ik. x e kA ] e" 



kA 
zk-r 



Z^lcA - 1 ) = - 



1 H [— ik x e kA ] e 



2 V 2A/7V 



2 V 2MiV 



'kA 



^k A 



X h [—ik x e kA ] e 



- ik-r 



2 V 2MiV 



V^kA + 1 



M;(k) = ( VkX 



1 / ft [ik x e kA ] e' 



-ak A 



^ kA + 1 > = - 



1 h [ikx e kA ] e l 



2 V 2MiV 



2 V 2AfiV ^/Z^ 
In Eq. one has, e.g., 

= E^T E e-^-^(6 a ,-e 7 + M(3™'M!(k)) (s« • M_(k) 

k A i^kA,-- 



-V^kA + 1- (211) 



Y— y 

kA D r/icA,-. 



e"^- /(kflT) « - £7 - ^k) (S« • M;(k) ) (3« • M+(k) 



After averaging over phonon populations, 



- E 



1 



e Ku k /(k B T) _ 



(212) 



(213) 



this becomes 



— y 



8MN 



kA 



„ (1) [k x e kA ] \ / (1) [k x e kA ] 



■'7a' 



(5 (e a / — e 7 + frw k ) n k 



„m [kxe kA ]\/ (i) [kxe kA ]\ 



(214) 
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This can be further simplified using 

[k x e ktl ] = ±fce kt2 (215) 
and, for the summation over the two transverse polarizations, 

£ (e kt -a)(e kt -b) = (a-b)- (k - a ^ k - b) (216) 

t=ti,i 2 

with a = b = e Xl and then e a . Averaging over the directions of the vector k: 

((k-a)(k.b)) = y(a-b) (217) 

one obtains 

k 



k 

= (SW • S«) £ ^ [<5 (uv 7 + Wk ) n k + S (av 7 - Wt ) (n k + 1)] . (218) 

k 

Now recalling Eq. (IMt one obtains 

(1) = irD* v fc 2 

a/3,a'/3' \2hMN ^ w k 
k 

Qia', 77 [<5 (% 7 + Wk) n-k + <5 (uv 7 - w k ) (n k + 1)] V/3 

-<W ^ Q^,77 ( W /3'7 + W k) ™k + 5 (Wp>y - Wfe) (n k + 1)] 

7 

/?)]}, (219) 



where 



is a dimensionless combination that characterizes the spin. Next, it is convenient to go over from summation to 
integration, 

1 v-^ f d 3 k . 
— \...=>v / 221 

where i>o is the unit-cell volume. Using vq/M = X/p and cjfc = Vtk one can introduce the characteristic frequency fit 
and the corresponding energy E t of the spin-phonon interaction 

n t =^y /4 , ^(p^ft 3 ) 174 . (222) 
As a result one obtains the characteristic relaxation rate 

vj £ 4 ( „„ _ „ 0) . ^ *M . r»M, (223, 



where 
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that enters the relaxation terms. In terms of T^^coq) and 



one obtains 



n,,, = 



z hu/(k B T) _ i 



R «L' > = -E^,77[ r(1) K Q 'K^+r (1) K',)(^+i) 

7 

-L^QSi 77 [rW(u^)n Uy0 , +r( 1 )(^ 7 ) (n u ^ + l) 

7 

+$43,0* p [ r(1) ("<*> a ) + !) + r(1) ("W) n^, + (« -» 0) 



(225) 



(226) 



Remember that here all r W (a;) with a; < are zero. Here Q^j q//3 , is defined by Eqs. (f220|) and (|202"|) . 

Within the secular approximation, all relaxation terms in the DME are defined by T aa i — R aa ,a'a'\ a i^ a , see Eqs. 

(f72|) . (f75|) . and l|77p. whereas for the direct processes considered here one has r^j = 0. For T^, from Eq. (|226l) one 
obtains 



r (1) , = 2 



1 'aa' 



_D 2 



r« ( Wa , a ) + 1) + r« (uw) ■ 







)- 





where we used 

/c« . bW ^ = / wW . o«* 

To compute the matrix elements of the spin operator components above, one uses 

1 i 
S± = S x ±iS y , S X = -(S-+S + ), S y = -(S--S+) 



S-l 



S-l 



S- = \m) l„ hm+1 (m + 1\ , S + = |m+l)i m+ i, 



m=-S 



m=-S 



m— — S 



where i m , m ' = y ^(5 + 1) — mm' . Thus one obtains 



(227) 
(228) 



5 2 = E |m) m(m| , (229) 



(a\S z \a'} 
(a\S x \a') 



(a\S y \a') 



y~^ (a|m)m(m|a') 



m=-S 



^ ((a|m)Z m , m+ i(m + l|a') + (a|m + l)l m+lt , n (m\a')) 



i 
2 



m=-S 
5-1 



(m|a')) 



(230) 



m— — S 



4- Raman processes 

Raman processes arise due to in Eq. (|196[) . as well as due to in the second order of the perturbation 
theory (details can be found in Ref. [I3|). The latter terms are nonessential at higher temperatures where Raman 
processes can become non-negligible since they contain a large thermal phonon frequency in the denominator. In 
Raman processes a phonon k is absorbed and a phonon q is emitted or vice versa. Processes with emission or 
absorption of two phonons make a small contribution and they will be ignored here. Thus the relevant phonon matrix 
elements are of the form 

(^kA - 1, fqA + 1 \8(j)^Slj> e | 2/kA, 
= (Vq\ + 1 ^qA) (^kA - 1 \5<f>£'\ ^a) + (^kA - 1 \S(/) ( \ (l>q\ + 1 \S<f>£'\ 

= M_,£,(k)M +i4 (q) + M .,;ki.\/..,Mq: = 2M^(k,q), (231) 
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where the matrix elements M(k) are defined by Eq. (|2 1 1 [) . Similarly to Eq. (|212[) . for one of the parts of R a p a ip> 
one obtains 



_Lr e -M(feT) 5 / +£ -e ,)v {2) v {2 \ 



b 

zuzu 



kA.qA b ^kXi^qX " 

X (s^M^k.q)) ( S ?icC^CC'(k,q)) 

J E <K E «' " £ 7 + - fiwk)s£J ift ,i (M* ir (k)M| >e (q) + M^(k)M^,(q)) 

kA,qA 

1 

7 q',CC' 2 

I ]T 5 ( Ea , - e 7 + fi^ q - 7^ k ) si 2 4 f/ M* !e (k)M; )e ,(q)H^ £r M_, c (k)M +iC (q), (232) 



x^ 2) , „,t (M_ x ,(k)M +x (q) + M_ iC (k)M+, c /(q)) 



ft 

kA,qA 



where 



n( 2 ) 

"q 7 ,^C 2 



2 ( S "7,??' + 5 «7,£'«) ( 233 ) 



and M are averaged over the phonon populations, j'kA — > "kA, see Eq. (|213[) . Substituting the expressions for M from 
Eq. (|211[) . one obtains 



T ( gj^]y ) E 5 ( £q ' ~ £ 7 + ~ ^k) ™k (flq + 1) 3 

^ ' kA,qA 

[k x e kA ] c [k x e kA ] c [q x e qA ] e , [q x e q x] f , 



(234) 



Now, using Eqs. (|215[) - (|217p , one can simplify this result to 

7T . \ , l1 n-(2) -(2) f 2 \ 2k2 1 2 x x 

(8MiV) 2 k f^ A 17 q 1 V q j 7Q ' CC W ^ k CJ q ?c £C 

= (l^iVf E — ^ + -q - "*) »k K + 1) • (235) 

Now recalling Eq. (|64[) one obtains 

(2) _ ttD 2 ^ fc 2 q 2 



R a lai = n 7 ny- (n n + 1) 

afi.affi (l2MN f j£ ^ k ^q ^ q ^ 

X 1 ~ E ( Wq '7 + W q _ w k) <W - E ^'kr/ (^'7 + U 1 ~~ Wk ) 

I 7 7 

+ Qal,a'P' [ S ( W W + ^q - ^k) + <5 {U3 aa t + UJq - Wk)] | , (236) 



where 



)(2) 

J a/3 ,ct' (3' 

a 1 



_Lv= (2) ^ (2) (237) 



similarly to Eq. (|220[) . Since Raman processes can become important only at high temperatures, one can drop spin 
transition frequencies in the energy 8- functions. Next, one can replace summation by integration with the help of Eq. 
(|22ip and introduce the characteristic Raman rate 

r(2) = 24M / ^ k " k K + 1} ' (238) 
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(2) 

Then R L , a , can be written in the form 

afj,a p 



7 

In the secular approximation one needs the rate 



- r( 2 ) (- 53 QW^Srp - <W 53 Qf Pm + agg,^ ] • < *w i 

\ 7 7 / 



r ( 2 ) d( 2 ) ! = 2r (2) o (2) = 2r(2) V ^ (2) ^ (2) 

act' aa.a'a' | . , , ^aa.a'o' £)2 / / "aa f ,^ /LJ a'a.^' 

_ fl W. _E!!lv|- (2) f (240) 



where from Eqs. ([205]) and ([233]) follows 



i 53 (e Q + e p - 2e 7 ) ((a |^| 7> <7 \S e | P) + (a \S e | 7) (7 \S £ \ /?)) 



7 



+ ^.9Ms (« l^| j8> + (a |%| /?)) - % 5~>Mb#«" ( a l 5 f"l $ ■ ( 241 ) 

Integration in Eq. (|238J) is limited by the Brillouin zone, so that Wk does not exceed some maximal value. We will 
use the Debye model in which the phonon spectrum continues in the same form up to the Debye frequency J7d that 
is the upper bound of integration. Thus Eq. (|238[) can be represented in the form 



24 2 7r 3 n Ef °\T 

where 

rv 

Aft, 

(e x - 1) 

For T <C Od the integration can be extended to infinity. Using Gq(oo) = 167r 6 /21, one obtains 



G n {y) = I dx — (243) 



r (2) - ^ (244) 



On the contrary, for T 3> 6d one can use Gq (y) = y /5 that yields 

r (2) = D2 ( fcB9p ) 5 ( fcsT ) 2 

2880tt 3 /L 



(245) 



The transition between these two regimes takes place at T/6 D ,i = l/y « [21/ (5 x 16tt 6 )] 1/5 w 0.2, i.e., much lower 
than the Debye temperature. For this reason the contribution of Raman processes is small in comparison to that of 
direct processes up to very high temperatures. Indeed, the contribution of Eq. (|244[) could become essential at high 
temperatures but long before it could happen the growth slows down to Eq. (|245|) . Then from Eqs. (|245p and (|223p 
one obtains the ratio 

r (2) s 1 (fc B e D ) 5 (fc B r) 

V rWcoth^a, 240^ 2 (hio ) 2 Ef ' 1 ' 

For Mni2 with 8d — 30 K, Et/ks — 150 K, and hajo/kB — D ~ 0.66 K (near the top of the barrier) the ratio rj = 1 
requires T = 2 x 10 4 K. 



D. The realistic phonon spectrum 



In the above derivations we have assumed that the crystal lattice possesses two degenerate transverse phonon modes 
that contribute to the spin-lattice relaxation. This is only the case for isotropic elastic bodies. In real crystals all 
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three acoustic phonon branches are different and they are neither fully longitudinal nor fully transverse. Fitting the 
heat capacity data for Mni2 to the extended Debye model of Ref. LUj has shown three acoustic modes with speeds 
Vi = 1541 m/s, v 2 — 2488 m/s, and v 3 — 3176 m/s. Since r' 1 ' oc v7 , only the mode with the lowest speed should be 
taken into account. This mode can be considered as approximately transverse. Thus we introduce the factor 1/2 in 
Eq. (|223p and in all other formulas for the spin-lattice rate due to direct processes. Similarly, the factor 1/4 must be 
introduced in Eq. (|238|) and in all subsequent formulas for Raman processes. The values of 0d and Et quoted below 
Eq. (|246l) correspond to v\ . 



E. Ground-state tunneling and relaxation 



1. The two-level model 



Consider the case in which H x and H y in Eq. (| 186[) are small, so that, in the absence of tunneling, the spin 
eigenstates \a) are basically \m) that are only weakly hybridized with the states in the same well. We will denote 
these states by 

s 

IVO = E <W'|"*">, (247) 

m"=-S 

where |c mm | = 1 and all other coefficients are small. Near tunneling resonances, these states are strongly hybridized 
with resonant states in the other well. Hybridization of the states \ip m ) and |Vw) can be taken into account in the 
framework of the two-state model 



tp } = e m t = to, to' 



<V) = (248) 



where A is the tunnel splitting of the levels to and m! that can be calculated from the exact spin Hamiltonian Hs 
or determined experimentally and ip is a phase. Since one can multiply the basis functions \ip m ) by arbitrary phase 
factors, we will set (p = it for convenience. This will result in a simpler form of the wave functions than in Ref. 0], 
whereas all physical results remain the same. Then the model above can be formulated as the pseudospin model 

fleff = -~o--A+-(e TO +e w /), (249) 

where components of <r are Pauli matrices, 

A = Ae x + We z , (250) 
is the effective field, and W is the energy bias or resonance detuning 

W = e m -e m ,, (251) 
defined by Eq. (|189p . We will need the direction angle 9 of A, 

W 

cosO = —= . (252) 

VW 2 + A 2 V ' 

The pseudospin acts on the states as 

O-x \lp m ) = IVV> , O-Jt \lp m ') = \lp m ) ■ (253) 

Of course, one also can calculate matrix elements of the physical spin S with respect to this basis. 

Eigenstates of H e g are the states polarized parallel and antiparallcl to A, and the eigenvalues are given by 

e a = \ (e m + e m , + a\fw 2 + A 2 ^j , a = ±. (254) 
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The transition frequency ujq between the levels is defined by 



e+ - e- = hu = y/W 2 + A 2 . (255) 



The eigenstates of H D g can be expanded over the natural basis as 



I") = \Xa) = ^ (C a IVO - aCU |^ m ,» , (256) 



where 



C Q = Vl + acos6». (257) 



One can directly check H e g \x a ) = s a \x a ) ■ The coefficients C a satisfy 1/C a — C_ Q / sm8. 

Let us prepare the system in the state \ip_ s ), the ground state in the left well that can be at resonance with 
the ground or excited state in the right well, \if> m /) with m! < S. Near the resonance these states hybridyze into 
V'i)* At low temperatures, all levels above i/*— s) are unpopulated and do not contribute to relaxation. The only 
relaxation processes are between \ip±) and the levels in the right well below |Vw) • Again, since the levels are only 
weakly hybridyzed inside the wells, here the dominant process is decay \ip m >) —* |Vw+i) ■ The inverse process can be 
neglected at low temperatures. In the case of the ground-state resonance, W-s) with \ipg) , this decay process is, of 
course, absent. The full description of both ground-ground and ground-excited resonances at low temperatures includes 
only two levels \x±) ; so that the effective DM is a 2x2 matrix. The DME in the general case of time-dependent spin 
Hamiltonians is Eq. (f6"3")) with non-adiabatic terms from Eq. (|55|) added, i.e., 

d P^fi = 1*7 ) p iP + PaiiX-y \kp)) ~ iu aP p aj3 + R a /3, a >0'P a >(3', (258) 



dt 



a'P' 



where all indices take the values ±. If \ip m i) is an excited state, in general one cannot use the secular approximation 
because loq can be comparable with the relaxation rate. 

Let us work out the non-adiabatic terms in Eq. (|258[) . Calculating the time derivative of C a in Eq. (|256[) , 

da= 2^^ COS ^~2^ Sin ^ = _aC -4' (259) 



one obtains 



\ Xa ) = i (C a \4_ s ) - aC_„ |</w>) = ~ |^_ s ) + aC a \^ n ,)) °- = -a | X _ a ) t. (260) 

Thus in Eq. (|258|) the scalar products are 

6 

(X a \X(i) = -a-{X- a \xp) = -a5- a ,p-. (261) 

and 

Q 

((Xa \Xj) Pyf3 + P al (Xj \Xp)) = ~ { a P-ot,fi + PPa-p) 2 ■ ( 262 ) 

7 

The density operator in the initial state typically is 

P(0) = \Tp- S )(i>s\, (263) 
so that the density matrix in the diagonal basis is given by 

P*p(0) = {a \p(0)\ (3) = (a\ V_ s ) (^_ s \ P) = \c a Cp, (264) 

where Eq. (|256|) with m = —S was used. In particular, 

P++(0) = I (I + cos 9), P __(O) = i(l-cos0) 

P+-(0) = p_ + (0) = isinfl. (265) 
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2. Ground- ground state resonance 



The results obtained above already allow to consider the dynamics at the ground-state resonance, m! — S. In this 
case the relaxation terms in Eq. (1258j) contain only T{uiq) lhq, so that the secular approximation is applicable. 
Dropping nonsecular terms in Eq. (|258|) one obtains 



d 



- {p+- + P-+) 2 + R++,++p++ + R++.--P- 

- [P— - P++) 2 _ tuJ aP+- + R+-,+-p + _, 



whereas 



Using Eq. (|226p . one obtains 



1-Ph 



(p + - 



-2Q++,__r(o;o) (n Uo + 1) = 
-2Q__, ++ r(w )n Wo = -r+_ = _ e -fi-o/(fe B T) r _ 



R —,- 

R—,++ = 2Q__ >++ r(w„) K„ + 1) = r_+ 



_ -feo/(teT), 



and 



#+-,+- = -<2++, — r(w ) (n Wo + 1) - Q__, ++ r(o;o)n tl , = -- (r_+ + r+_) 



etc. Thus Eq. ()266|) takes the form 



(266) 



(267) 



(268) 



(269) 



P++ = -(p+-+P-+)2- T -+P++ +T +-P- 



P+- = (P++-P-) 



This system of equations can be rewritten as 



iuj + ^ ( r -+ + r +-) 



P++ = ~ e - r - P++) 

p+- = e(p ++ -i/2)-(tu + r/2) P+ _, 

where 

r = r_+ + r+_ 

is the total relaxation rate between the ± states and 

P e -hu /(k B T) 



eq 
P++ = 



(270) 



r_+ + r+_ i + e -^o/(fcBT) 



(271) 



(272) 



(273) 



is the equilibrium population of the upper level. 



3. Dynamics of the ground-ground state resonance via effective classical spin 

The DME for the ground-ground state resonance, Eq. (|270p . can be conveniently formulated in terms of the averages 
the pseudospin <x with the density operator. Using Eq. ((¥]) one can write 



a = (*) = XI P c0 (X/3 Xa) 
a/3 



(274) 
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Directing the axis z' along the total field A, one has 

e> = (ct_ + 0+) e x > + i (<r_ - <t+) e. y - + a z >e z > 

= (\x+) (X-\ + \X-) (X+|) «V + i (x-| - \X-) (x+\) e V + (I*-) ~ <V- (275) 



Then one obtains 



OV = p + -(X-\vx>\x+) + P-+(x+\vx'\X-) = P-+ + P+-=2R£P-+ 

<ty = P+-{X-\v y >\X+) + P-+{x+\Vy>\X-) = i (p-+ ~ P+-) = 2Imp_ + 

oy = P++(x+\vz'\x+) +P — (x~ |ov|x_) =p — = 1 ~ 2 P++- ( 276 ) 

Now Eq. (|270p can be transformed as 

ov = /0_+ + />+_ = {p ++ - P ) + iuj (p_ + - p + _) + R x , = -9<r z > + u) (7y> + R x > 

&yi = i {p_ + - p + _) = i (iuj a p_ + + iuj a p + _) + R y f = -ujq<t x > + Ry< 

g z , = ( P+ _ + p_ + ) 9 + R z ,= ho x , + R z , (277) 

<T = [ffx(wo + fi)]+R, a; =cj e Z ', fl — 9e y i. (278) 

This is a Larmor equation for the classical vector er in the frame rotating with frequency Jl due to the time dependence 
of the spin Hamiltonian. The relaxation vector is given by 

R = -i(r_+ + r+_)(7^e 1Ef -i(r_ + + r + _)(7 1 , f e v /-2(-r_ + p ++ + r + _ / >__)e,/ 

= ~ (r_+ + T+_) a x ,e x ,~ (r_+ + r+_) a y ,e y , - [-T-+ (1 - a z ,) + T+_ (1 + o z <)\ e z , 

= -\ (r_+ + r+_) a x ,e x ,-\ (r_+ + r+_) a y ,e y , - [(r_ + + r + _) g z , - (r_+ - r+_)] e z , (279) 



or 



or 



R = -- cr — w - T cr cq , (280) 

2 V w n / V 



where T is given by Eq. (|272|) and 



r_ + -r + _ fiw 



= tanh— ^ (281) 



r_+ + r + _ 2fc B r 



is the equilibrium spin polarization. 

Eq. (|278p can describe, in particular, the Landau-Zener effect of transition between the energy levels of a two-level 
system as the energy bias W is swept though the resonance, W = 0. It is remarcable that this essentially quantum 
phenomenon can be described by a classical language. Studying the dynamical behavior of the classical spin a helps 
a lot understanding the LZ effect. In particular, if the sweep is slow, the non-adiabatic term ft in Eq. (|278[) is small 
and the spin remains nearly collinear to ujq at all times. In the adiabatic frame u)q — const and the spin vector in the 
laboratory frame (i.e., in the natural basis) can be obtained by a rotational transformation. One also can rewrite Eq. 
(|278jl in the laboratory frame by simply dropping fl and making u>q time dependent according to Eq. (|250|) . = A. 
Because of the vector form of the equation of motion, its transformation between the frames is easy, in contrast to the 
transformation of the DME between the natural basis and the diagonal (adiabatic) basis described in Sec. Ill Gl As 
the vector A is rotating, cr is lagging behind it, depending on the sweep rate. For a slow sweep rate it nearly follows 
the direction of A while for a fast sweep it nearly remains in the initial state. 



4- Coherence in the ground-ground state resonance 



Let us consider the time evolution of a two-level system in the case of time-independent spin Hamiltonian, 9 = 0. 
This can be done using either Eq. (|271[) or Eq. (|278[) . In particular, the solution of Eq. (|271|) with the initial 
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conditions of Eq. (]265[) is 



P++ (t) = p c + \ + [p ++ (o)-p c ;\ 

P+ _{t) = p + _(0)e-< to "+r/a>. = l 



-VI 



P++ 



W 



VW 2 + A 2 



eq 



-(tw +r/2)t 



(282) 

2 V^ 2 + A 2 V ; 

Transformation back to the natural basis is done using Eq. (|86[) . The probability to remain in the initial state V 7 — s) 



is 



a/3 
1 

2 
1 
2 
1 



IF 



Viy 2 + 


AV 


w 


A 2 ) 


Vw 2 + 






A 2 ) 


V^ 2 + 




A 2 





p++(t) 



a/3 
1 



1 - 



w 



w 



y/W 2 + A 2 



) [1 -^ (i)1 + vra Rep+ - (i) 



or 



P-s.-si 1 ) 



2W 2 + A 2 
W 2 



Vw 2 



VW 2 + A 2 
"1(1 



Re P+ _(t) 



w 



VW 2 + A 2 



eq 



,-rt 



(283) 



-r< 



VF 2 + A 2 
W 



w 2 



-W cos (w t) + 1 



-rt 



Vw 72 



-rt\ 



(284) 



VM^ 2 + A 2 " ++ 

The time evolution described by this formula is a combination of relaxation with rate T and oscillations of frequency 
wo damped with rate T/2. It should be noted, however, that the ground-state tunnel splitting A is typically very 
small, so that it is very difficult to experimentally realize \W\ < A to see coherent oscillations of the spin between the 
two states. 



5. Relaxation rate between two tunnel-split states 

The relaxation rate V |_ for the ground-state doublet can be found analytically [8]. In particular, for the uniaxial 

model in the presence of transverse field along the x axis, with the help of the high-order perturbation theory one 
obtains 

A ml — m 



(<>P-\S x \i> + ) 



VW 2 + A 2 2 
A W ml — m 



9Pb h x VW 2 + A 2 



Then from Eq. (|202|) in components, 

E-+,x = -itiua ("0- \S x \ip+) - \S y \ tp + )g^ B H z + (tp_ \S z \ip + ) gp, B H y 
B_+,„ = -iftu> (ip- \S y \ip + ) - ("0- \Sz\ip+) 9P-b h x + (ip- \S X \ ip + )gfi B H z 
S_ +iZ = -iriuj Q (ijj_\S z \ilj + ) ~ (ijj_\S x \ijj + ) gp, B H y + (ip_\S y \ip + ) gp, B H x , (286) 



with H y = one obtains 



^- + .2 



A ml - m 

(W - gp B H z 



gp B H x 2 



- A m '- m rr ( -i A 2 W^-g^g,) ^ 

~ V^ + A 2 2 ^ I 1 - {9 , bHx) 2 ) • < 287 > 
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Taking into account A <C g[i B H x and using W = (m' — rn)g[i B H z yields 



(m'-m- \f^l + ^ Hx }l (l-(m f - m)(m' - m - 1)^1) 



Now, from Eqs. (J255J), (l2"2"0j) , and (|2"2"5f one obtains 



(288) 



r_+ = 2Q__ )++ I> ) = |S__ 



m' — m\ 2 A 2 l<jo (g^ B H x 



12nE? 



W 2 + A 2 
(gH B H x y 



H ( H 
(ml - m - l) 2 ^- + ( 1 — fm' - m)(m' - m - 1)— | 
Hi V Hi 



(289) 



Dropping again the small A 2 term in the square brackets one obtains 



r_ + = 



ml - m\ 2 A 2 uj (g/j, B H x Y 



12tt Ef 



(290) 



where 



Q = £ 4 + (i-£ 2 ) 2 = i-2£ 2 + 2£ 4 

TT 

£ = ^/ (m' — ml) (ml — m — 1) — -. 



(291) 



In particular, for the ground-state resonance, m — —S, ml = S, one has £ = y/2S(2S — l)H z /H x . One can see that Q 

begins to deviate from 1 only for SH Z > H x that requires W ^> A. For such a strong bias, the dependence of F |_ on 

the bias has the form 



T_+ ex /(£) = CQ - £ (1 - 2£ 2 + 2£ 4 ) . 



(292) 



The function /(£) monotonically increases but its slope initially decreases from 1, attaining the minimal value /'(£ ) = 
1/10 at £ = y/3/10 ~ 0.5477. At this point /(f ) = ^3/10(29/50) ~ 0.3177. After that the slope of /(£) begins 
strongly to increase, so that T |_ ~ £ 5 ~ H z . 



6. Ground- excited state resonance 

Let us consider the tunneling resonance between jV'-s) ano - IVw) with ml < S. The hybridyzed states \tp±) 
decay, predominantly, into |Vw+i) • The corresponding relaxation rate is large, so that we will neglect the rate 
of transitions between |V'±) that is small near the resonance. At low temperatures one can neglect the energy-up 
processes \ip m '+i) ~ * \' t P±) ■ Thus, as in the case of the ground-ground state resonance, it is sufficient to consider the 
2x2 DME for the states \ip±) ■ I n this case, however, the normalization of the effective DM of the two-level system is 
not conserved because of the decay to the lower states. 

The DME for the states \ip±) can be obtained from the general formalism, as above, but this way is lengthy and 
the final results can be written without any calculations. In fact, the system obeys the damped Schrodinger equation 
in the natural basis that is simpler than the DME and has the form 



d 


i W 


i A 


dt C - s = 


~2~h c ~ 


s ~ 2~h Cm 


d 


(i W 


v 


di Cm ' ~ 


V2~T " 


2" 1 - m' + l.m' 



(293) 

where the level ml is damped and the level — S is undamped. The decay rate between the adjacent m-states in the 
generic MM model is given by Eq. (A9) of Ref. Q that can be rewritten as 

_ (2m+l) 2 l 2 n+hm (D/h) 2 \u m+1 , m \ 3 
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Here l m +i,m is defined below Eq. (I229|) and il t is denned by Eq. (j222[) . In the case m — S — 1 and m + 1 = S, Eq. 
(|294j) simplifies to the elegant form 

S 2 u% i c 

where huis-i,s = — 1) D. 

The DME can be obtained from Eq. (|293| by setting 

P-s-s = c -sc*^ s , p m , m , = c m ,c* m ,, p_ w = c_ s c^, (296) 

and calculating time derivatives. It has the form 



2h [P - s >* 



df rm ,m 



Pm'-S) 

i A 
' m ' ~ 2~ft 



{P~S,m> ~ Pm'-s) 



,W 1 \ iA, . 

Q W+l.m' I P-S,m> + 2^ \P-S-S ~~ Pm'.m') 



(297) 



that coincides with the results of Ref. (In the latter the precession goes in the wrong direction, however). It 
should be stressed once more that this tunneling DME is non-secular. 

Of course Eq. (|293p is easier to solve than Eq. (|297[) . We search for the solution of Eq. (|293[) in the form e~ At . 
Eigenvalues A of Eq. (|293[) satisfy the equation 



or, with the use of Eq. l|255j) . 



that yields 



i_ A 

2 H 



2 



_ A 

.2 h 



= 



-A + 



i W 
2~h 



i_W_ 1 \ 1 / A 

2 ft + 2 raW + 4fi 



a — ~r m '-i-i im 'A 



iW 
4~h 



■ m' + l,m' 



(298) 



(299) 



A± 



1/1 /l .W 

2 [ 2 r,n ' +1 ' m ' 1,1 V 4 rm '+ 1 > m ' ~ w o ~ l y r ™'+i, 



2 V 2 

The solution of Eq. (|2"9"3")) has the form 



/ + \ _ (+) -A + t , )„-A_i 



The eigenvectors (a^,^) follow from Eq. (f293|> : 



A+a 



(±) 



±a_s 



2~fi~' 



(±) 



A 



(±) 



2 ft° m ' 



(300) 



(301) 



(302) 



i.e. 



,(±>_ h(w_ 
A I h 



(303) 
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Thus Eq. (|30ip can rewritten as 



c_ s (t) = aL + 2e- A +* + aL"ie- A -* 



c m '(t) 



(+) W/h+2i\ + _ x { _ ) W/h+2iX- _ A _ t 

-a_s 7-77 e — a c T-77 e 

6 A/ft ~ 6 A/ft 



From the initial conditions c_s(0) = 1 and c m '(0) = one obtains 

(W/h+2i\ + ) + a ( ~<! (W/h+2i\_) = 
(+) _ 1 _ W/H+2iX- _ iW/(2h) - A_ 



that yields 



1 w/h+2ix+ 2i(A_-A+) A+-A_ 



( _) _ W/(2fi) - A -1 



and, from Eq. l|50g|) . 



A- - A, 



(+) _ 9 , ft W(2») - A+) (»W/(2fi) - A_) 



Finally one obtains 



c-s(t) 
c m > (t) 



,(-) 



1 



A 

,(+) 



i W 

2T~ X - 



A+ - A_ 

2i h fiW 



\+ - A_ 



,-A+t 



2I1 H 



,-A_t 



A+ - A_ A V 2 ft 
Now the components of the density matrix are given by 



2l 



A_ e- A + 



1 



IA, - A I 



i W 

2T~ X - 

i W 



IA+ - A_ 



2 Re 



2 ft 
2~ft ~ A - 
2l ~ 



,-A+t 



___ _ a* I e" A +* - [ --— - AI ] e 



i W 

2T _Ah 



A+t e -A_n 



,-A-t 



e -2Re(A + )t + 



i W 



2 ft 
i W 

2T~ X - 



,-A* t 



-2Re(A_)t 



A * e -(A ++ A*) t 

2ft + 



and 



A 



2 a\±W_ 

* I 2 ft 



A I 



2T _ A +l 



IA+-A- 



-2Re(A + )i 



e -2Re(A_)t_ 2Re ^-(A + +Al) t -J 



(304) 



(305) 



(306) 



(307) 



(308) 



(309) 



(310) 



One can see that the relaxation is described by three relaxation rates, 2Re(A±) and Re(A+ + A_). In addition, there 
are oscillations with the frequency Im(A + — A_) corresponding to quantum transitions \—S) f± |m') . 
Spin polarization in our low-temperature tunneling process is given by 



(S z ) = -Sp_ s 



s 



S ' / , ""from - 
m—m' 



(311) 



As the states with m = m' + 1, . . . , S — 1 decay faster than \m') , their contribution can be neglected. Then one can 
write 

(S z ) = -Sp_ s _ s + m'p m , m , + Sp ss 

= ~SP-S,-S + m 'Pm>m> +S(1- P-S-S - Pm'm') 

= S-2Sp_ s _ s - (S-m')p m , m >. (312) 
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As both p_ s _ s and p m i m i go asymptotically to zero, the magnetization change in the process is (S z ) — {S z ) = —25, 
and 

(S z ) t - (Sz)^ = -2Sp_ s _ s - (S - m')p m , m ,. (313) 
Thus the integral relaxation time is given by 

n,t = {Sz)q _ {Sz)oo = J o * [P-8,-sW + -2g-/W«'(*)J • (314) 

For a small bias, say, m' — S — 1 and large spin, the contribution of the second term is small. Using Eq. Q309|) . one 
obtains 

iff \ I 2 iff \ I 2 /iff \ W iff 



dtp s s (t) = g j ^ ^ J + ^ ' -2Re-^ " j 1 * r ' ^ [> (315) 

I A . A I 2 1 2Re(A+) 2Re(A_) A+ + A! ' 



and 



S\ 2 4|if -A_|'||f -A + |" f I 1.2 



Mathematica-aided simplification yields 



and 



Then from Eq. (I314|) one obtains 



AW 4- A 2 4- r 2 

dtP-S.-SW = f T2 ( 317 ) 



dt Pm ,„At) = . (318) 

i m' + l,m' 



_ 4W + A 2 + r 2 m , +1 , m , 5 - m' 1 4W + (l + A 2 + r 2 „, +1 , m , 

Tint — r A2 oq v ~ r A 2 ' v^^^/ 

The corresponding rate can be written as 

In comparison to the overdamped result of Ref. [12j , this formula contains an additional term ~ A 2 in the denominator. 
In the underdamped case A 3> ftr m '+i,m' at resonance W = 0, the rate is given by 

Tint = (l + (321) 

that is of the order of thermal activation rate at high temperatures. This means that the barrier is completely cut at 
resonance and the relaxation rate does not significantly increase with temperature. Condition A 3> HT m i + i tm i makes 
the secular approximation valid, unlike the overdamped case A < HT m i+i jm / , where the secular approximation leads 
to unphysically large relaxation rates at resonance. 



F. Numerical implementation and illustrations 



In this section some representative results of the numerical solution of the DME for molecular magnets are shown. 
The code implemented in Wolfram Mathematica is available from the author. It is based on the diagonalization of 
the dynamical matrix $ of the DME for the free-evolution problem, as described in Sec. IIII Al Linear dynamical 
susceptibility, Sec. IIII C[ can also be obtained on this way. As the tunnel splitting A can be very small, especially in 
Mni2 (S = 10) at small transverse fields, matrix algebra requires using a high custom precision, typically 100 significant 
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FIG. 1: Zero-temperature escape rate from the metastable state vs the bias field in the generic model of MM. 




FIG. 2: Time evolution of (S z ) at over- and underdamped tunneling resonances in previous figure. 



figures. For N14 (5 = 4) standard precision is sufficient. Custom precision makes computation slower. Still, the DME 
in the secular approximation involving the (2S + 1) x (25 + 1) dynamical matrix solves fast enough on a standard 
PC. Solution of the full non-secular DME is very slow for S — 10 because of the big size (25 + iy x (25+ If of the 
dynamical matrix. It is very important to use the semi-secular DME that is no less accurate than the full non-secular 
DME but has the dynamical matrix of the size (65 + 1) x (6S + 1) . As a result, the solution, although slower than 
that of the secular DME, is still realistically fast. The difference between the secular and semi-secular versions of the 
method is confined to the close vicinity of the overdamped tunneling resonances, A < hT m i + i yTn i , while everywhere 
else the numerical results are the same. 

Numerical solution shows that the dynamical matrix $ of the non-secular DME has exactly 25+1 real eigenvalues 
out of the total (25 + l) 2 eigenvalues. One of real eigenvalues is zero and corresponds to the thermal equilibrium. 
Complex eigenvalues occur in complex conjugate pairs. This behavior is similar to that of the secular DME and 
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FIG. 3: Escape rate vs bias field in the generic model of MM at different temperatures. 



requires explanation. At low temperatures in the regime of thermal activation and weak tunneling (overdamped 
resonances) there is one nonzero real eigenvalue that is much smaller than all other real eigenvalues and real parts of 
complex eigenvalues. In the case of underdamped resonances, there are three eigenvalues, one real and two complex, 
that describe the slow dynamics. 

Fig. Q] shows the zero-temperature escape rate vs the bias field H z in the generic MM model with B = in Eq. 
(|188p . The striking feature is the spin tunneling at resonance fields that leads to the increase of the escape rate by 
many orders of magnitude. Most of the points have been obtained from the secular DME, the points at resonances 
and between them have been obtained from the semi-secular DME, and the analytical result of Eq. (|320p is drawn in 
the vicinity of resonances. Near the zero-field resonance, Eq. (I290|) is shown. The characteristic "shoulder" described 
by this equation is well reproduced by the numerical result. As mentioned above, the secular approximation can yield 
unphysically high escape rates at resonances but the resonances are narrow and there are no secular points that hit 
them. Resonances with k = 1,2, and 3 are overdamped and can be approximately described by the DME solution based 
on effective resistances method of Ref. [l2| • As explained in Ref. [12] , the barrier at resonances is lower than the full 
classical barrier, and its height corresponds to the so-called blocking level for which A mm / ~ h(Pm'+i,m' + r m _i -m ) ■ 
Resonances with k > 4 in Fig. Q] are underdamped, so that the blocking level is the ground state and the barrier is 
reduced to zero. Accordingly, the escape rate at these resonances is of order T ~ 3xl0 6 s _1 and it coincides with the 
spin-phonon rate between the adjacent levels that is the highest possible rate achievable off resonance at temperatures 
exceeding the energy barrier. 

Fig. [5] shows the time dependence of spin polarization (S z ) at different resonances in Fig. [TJ The relaxation is 
exponential for the overdamped resonances, as well as off resonance (not shown). To the contrast, at underdamped 
resonances with k > 5 there are damped oscillations described by three different relaxation rates in Eq. (|309[) . In 
the case of exponential relaxation, it is sufficient to identify the escape rate with the smallest real eigenvalue of the 
dynamical matrix. In the case of underdamped resonances there are three slow eigenvalues, and obtaining the correct 
value of the escape rate requires using the integral relaxation time. 

Escape rate vs the bias field at different temperatures in the generic model is shown in Fig. [3] All data were obtained 
from the numerical solution of the semi-secular DME. The anisotropy value D — 0.66 K has been chosen to fit the 
barrier height in Mni2 (see below). As expected, the escape rate increases with temperature, faster off resonance than 
on resonance. One can see (especially clear for T = 2 K and k = 1) that at nonzero temperatures the tunneling peak 
may consist of several peaks of different width on the top of each other pjj]. Broad peaks correspond to tunneling 
at high energy with a large splitting A while narrow peaks correspond to tunneling via a low-lying resonant pair of 
levels with small A. At zero temperatures the zero-bias tunneling peak is very narrow because of the anomalously 
small damping of the ground-state levels and it is not seen in the plot. However, for nonzero temperatures this peak 
becomes broad as all other peaks because of tunneling via excited levels that are regularly damped via decay to the 



44 




FIG. 4: Temperature dependence of the escape rate in the generic model on and off resonance. 




FIG. 5: Escape rate vs transverse field in the generic model of MM at different temperatures, on and off resonance. 



lower-lying levels. 

Arrhenius plot in Fig. [4] shows transition between the thermal-activation and ground-state-tunneling regimes on 
and off resonance for different transverse fields parametrized by h x = gfi B H x / (2SD) . For small transverse fields, the 
resonances are overdamped and ground-state tunneling is small. In this case the activation part of the plot is nearly 
a straight line with the slope corresponding to a particular effective barrier. The transition to the horizontal line 
describing the ground-state tunneling has little rounding. This is the so-called first-order transition between thermal 
activation and ground-state tunneling, the two competing channels [HjElj]- Fot h x = 0.2 at resonance, the activation 
part of the plot is noticeably curved. This is a manifestation of the second-order transition in which the dominating 
tunneling level gradually moves down with lowering temperature, effectively decreasing the barrier height and the 
slope of the curve. For h x = 0.3 at resonance, the barrier is reduced to nearly zero and the ground-state tunneling is 
very strong. 
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FIG. 7: Temperature dependence of the escape rate in Mni2. 



Fig. [5] shows the dependence of the escape rate on the transverse field h x at different temperatures on and off 
resonance. On resonance at nonzero temperatures, there are characteristic steps arising as a result of moving the 
blocking level up or down the energy [13]. This phenomenon can be seen in Fig. 5 of Ref. [13], obtained by the 
effective resistances method. One can see that on resonance, k = 1, the barrier goes to zero with increasing h x , so 
that above some critical value of h x curves corresponding to different temperatures merge at the level of the highest 
possible rate. At these transverse fields the barrier off resonance still exists since the curves corresponding to different 
temperatures merge at higher values of h x . 

Next figures show the numerical results for Mni2. Because of the quartic uniaxial anisotropy B, tunneling peaks in 
Fig. [6] are split, as explained in the comment after Eq. (|T91j) . The rightmost big peaks correspond to the ground-state 
tunneling, and smaller peaks to the left of them, seen at nonzero temperatures, are due to tunneling via excited states. 
Graphed results of earlier calculations of this kind for Mni 2 can be found in Refs. [3, [HI]. In comparison to the 
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results for the generic model with the same barrier 66 K in Fig. [4j Mni2 shows ground-state tunneling up to higher 
temperatures. 

Temperature dependences of the escape rate in Fig. [7]are different for different bias fields. If for a given H z there is 
a tunneling resonance at some energy between the top of the classical barrier and the ground-state, thermally activated 
tunneling via this resonant pair is a relaxation channel competing with the two channels considered above. As a result, 
there are two different slopes in the Arrhenius part of the plot, such as for k = 1.006. This value of k corresponds to 
the high blue peak in Fig. [6] that disappears at T = 0. 



G. Discussion 



Existing work on molecular magnets using the density matrix equation can be split up into two groups: (i) using 
the natural or m-basis and (ii) using the diagonal basis. In all known cases the DME is reduced to the system of rate 
equations for the diagonal DM elements, the level populations. Using the natural basis is justified if the terms in the 
spin Hamiltonian that are non-commuting with S z are a small perturbation. However, even a small non-commiting 
perturbation can severely distort the levels near the top of the barrier that are mostly inportant in thermal activation. 
On the other hand, tunneling via robuster low-lying levels at low temperatures can be well described perturbatively 
in the m-basis. 

In Ref. [l^ | the thermal activation rate of a generic MM was calculated in the m-basis in the absence of tunneling 
via the integral relaxation time. Tunneling has been taken into account in Ref. fl2| by adiabatically eliminating fast 
nondiagonal DM elements that amounts to using the high-order perturbation theory in calculating tunnel splittings 



17 1 . The resulting system of rate equations with resonance tunneling was solved by the method of effective resistances 
12| | using the idea of the solution of the Fokker-Planck equation at low temperatures in the classical case. Later the 



system of rate equations in the m-basis was employed in Refs. [ill EE GjI 

In particular, Ref. fl5| re peats the steps of Ref. [l2j using the realistic model of Mni2 with B ^ in Eq. (|187|) . 
A new element of Ref. [15| is the erroneous consideration of spin-phonon interactions leading to the spin-phonon 
coupling of the type D i y S\ + S 2 .) (e xx — e yy ), e aa being components of the deformation tensor, because of tilting the 
easy axis by transverse phonons at second order in small tilting angle Sip. This leads to nonexistent direct processes 
with changing m by 2. In fact, as we have seen above, second-order terms in S<p, Eq. (|198[) . give rise to Raman 
processes rather than to direct processes. The error made in Ref. [151 ], neglection of a part of Sip 2 terms that cancel 
the result, has been explained in Ref. [2(|. Nevertheless, the appeal of Am = 2 direct processes has been remaining 
strong, so that the relevance of Eq. (A12) of Ref. [15j for explanation of experiments on molecular magnets is still 
disputable. The recent examples are experimental works on Fes, Refs. [2l], [22[. Whereas in Ref. [2l[ Eq. (A12) of 
Ref. [15] is used with success, Ref. [22j states that direct processes with Am = 2 do not fit the data. On the other 
hand, for Fes these processes were shown to arise from rotations around the easy axis, the corresponding coupling 
constant being the transverse anisotropy E, see Eq. (B5) of Ref. 

Moving to the universal form of the DME proposed here can help to end the confusion about what to include into 
the relaxation terms, since the latter automatically follow from the spin Hamiltonian and there is no freedom to make 
a mistake. Moreover, in the diagonal basis used in the universal DME there are phonon-induced transitions between 
all exact energy levels, not only between the nearest or second-nearest neighbors. An example is the relaxation rate 
between the levels of the ground-state doublet, Eq. (|290p and Eq. (69) of Ref. Q. 

Among the works using the DME in the natural basis, in the secular approximation, are Refs. [3,[23]. The authors 
say that the advantage of this method is that spin tunneling is absorbed in the exact basis states. This is overall 
true, although the full dynamical description requires taking into account the decoupled nondiagonal DM elements, 
in addition to the system of rate equations that was used. Also in the case of weak tunneling (overdamped tunneling 
resonances) the secular approximation fails and results in unphysically large escape rates at resonances. This was 
explained at the beginning of Sec. IIV C[ as well as in the comments below Eq. (19) of Ref. [HI and above Eq. 
(2) of Ref. [22I ]. Certainly something is missing if spin tunneling is automatically incorporated into the exact basis 
states but, in spite of it, one cannot approach tunneling resonances. The solution is to use the non-secular or better 
semi-secular DME that takes into account the dynamical coupling between the diagonal and slow nondiagonal DM 
elements and is thus valid everywhere. It should be stressed that the system of rate equations with tunneling in the 
m-basis is essentially non-secular and for this reason it does not fail at resonances. 
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